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ABSTRACT 

Six  degree- of- freedom,  rigid  body  equations  of  motion  are  described  suitable  for 
modeling  the  dynamic  characteristics  of  multistaged,  free- flight,  ballistic  rockets  such  as 
the  DRES  developed  aerial  targets  CRV7/BATS  and  ROBOT-9.  These  equations  of 
motion  form  the  core  of  a  FORTRAN  simulation  software  package  called  BALSIM. 
This  package  allows  for  modeling  of  vehicle  thrust  and  structural  asymmetries,  time- 
varying  mass  and  inertia  characteristics,  variable  wind  conditions,  nonstandard 
atmospheric  conditions,  stage  failures,  and  different  rocket  motor  types.  The  BALSIM 
package  has  been  written  in  IBM  FORTRAN  IV  and  has  been  tested  on  the  IBM  3033 
computer  with  the  H-extended  compiler.  It  is  currently  being  adapted  for  use  with  the 
VAX1 1/780  and  Honeywell  DPS-8/70C  computers. 
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vehicles) 
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Launch  rail  guide  length  (see  Figure  5) 

Thrust  vector 

Atmospheric  temperature 

Thrust  of  the  i-th  rocket  motor 
Components  of  V  in  F« 

Components  of  V£  in  F* 

Components  of  W  in  F* 

Airspeed  vector 

Velocity  vector  of  the  vehicle  centre-of-mass  with  respect  to  F, 
Magnitude  of  V 
Vuj,2  +  w  b1 

Wind  velocity  vector  with  respect  to  F / 

Components  of  W  in  F, 


(ix) 


UNCLA.  IHEr 


UNCLASSIFIED 


w //w 


(*ab,  y 


Ag9 


Z  Aa) 


(Xr„,  Y r„  Ztb) 

(X«c»  me 9  Zac) 

y 

(xCf,  y,*,  zCf) 


x  (L«i  y«m>  Z,m) 

(*r.  yr>  ZT) 

(X/,  y/f  Z,) 

[(X  (y«.)i,  (z*f«)i] 


(*mi  y«,  zM) 

[(x«t)<,  (ym)ii  (Zm)i] 


Pseudo  fin  airspeed  component  normal  to  its  chord  plane,  see 
equation  (2.4,11) 
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Components  of  Rr  in  FT  [see  equation  (2.10,3)] 
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Coordinates  of  the  centre- of- mass  of  the  empty  motor  case  of 
the  i-th  rocket  motor 

Coordinates  of  the  payload  centre-of-mass  in  F* 
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Angle  of  attack  of  vehicle  [see  equation  (2.4,6a)] 

Angle  of  attack  of  pseudo  fin’s  chord  plane  [see 
equation  (2.4,12b)] 
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Euler  bank  angle  for  FB 

Cylindrical  coordinate  of  pseudo  fin  (see  Figure  4) 

Euler  azimuth  angle  for  Ffl 

Angular  velocity  vector  of  Fa  with  respect  to  F, 
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1.  INTRODUCTION 

In  1979  DRES  began  development  work  towards  modifying  the  U.S.  army  BATS 
target  for  use  with  the  higher  specific  impulse  CRV7  rocket  motors  under  the  auspices  of 
a  joint  US-Canadian  TTCP  agreement.  This  original  work  has  lead  to  a  number  of 
DRES  initiated  activities  including  the  modification  of  the  target  to  permit  multistaging 
and  the  development  of  an  all  CRV7,  multistaged  vehicle  referred  to  as  ROBOT-9 
(Figure  1).  As  well,  support  equipment  has  been  developed  for  target  operation  in 
moderately  heavy  seas,  i.e.  in  conditions  typical  of  Canadian  coastal  waters  (the  ROBOT 
System  development  history  is  summarized  in  more  detail  in  Reference  1). 

To  support  the  development  of  these  free-flight  targets,  computer  simulation 
programs  were  required  that  predicted  the  dynamic  characteristics  of  the  vehicles.  Of 
particular  importance  were  accurate  predictions  of  basic  performance  parameters  (e.g. 
range  and  flight-time),  wind  effects,  effect  of  nonstandard  atmospheric  conditions  and 
the  dynamic  effects  of  launching  from  a  moving  ship  on  a  finite,  nonzero  length 
launcher. 
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No  existing  software  was  available  to  DRES  that  performed  all  of  these  tasks 
while  conveniently  permitting  some  configuration  variation.  As  a  result,  in  the  period 
September,  1980  to  December,  1981,  a  six  degree-of-freedom  simulation  package  was 
written,  debugged,  and  tested  at  DRES.  This  package  was  applied  to  the  evaluation  of 
the  CRV7/BATS  and  ROBOT-9  performance  characteristics  and  safety-envelopes.  It 
was  coded  in  IBM  Fortran  and  has  been  used  on  the  IBM  3033  computer  with  the 
H-extended  compiler.  The  package  is  currently  being  installed  on  a  VAX  11/780 
computer  for  use  with  a  FORTRAN  77  compiler,  and  will  be  adapted  for  use  with  the 
Honeywell  DPS-8/70C  computer.  The  package  has  been  designated  BALSIM 

It  is  the  intent  of  this  report  to  provide  documentation  of  BALSIM  in  sufficient 
detail  to  permit  users  familiar  with  FORTRAN  to  run  the  program.  Chapter  2  develops 
the  dynamic  model  and  summarizes  its  limitations.  Chapter  3  describes  the  BALSIM 
package  in  general  terms.  Finally,  Appendix  1  is  intended  as  an  essentially  self- 
contained  userbook  for  the  package,  and  includes  a  listing  of  all  program  modules. 

2.  DYNAMIC  MODEL 

This  section  summarizes  the  key  features  of  the  dynamic  model  programmed  into 
the  BALSIM  package.  The  basic  six  degree-of-freedom  equations  are  derived  in  the 
following  sections. 

2.1  Fundamental  Assumptions 

Several  overall  simplifying  assumptions  have  been  made  in  the  derivation  of  the 
equations  of  motion.  They  are  valid  for  ballistic  rocket  vehicles  that  have  rigid  structures 
and  relatively  short  ranges,  i.e.  less  than  100  km  (50  nm). 

The  assumptions  are  as  follows: 

1.  The  Earth  is  flat  and  any  Earth-fixed  reference  frame  is  inertial. 

2.  The  vehicle  is  a  rigid  body. 

3.  There  are  no  control  surfaces. 

Assumption  3  may  be  readily  relaxed  by  adding  the  appropriate  control  terms  into 
the  equations  of  motion. 


UNCLASSIFIED 


UNCLASSIFIED 


n 


2.2  Reference  Frames,  Rotation  Matrices  and  Angular  Velocities 

In  the  general  case  both  an  Earth-fixed  (say  F£)  and  an  inertial  reference  frame 
(say  F,)  must  be  defined.  Because  of  the  first  simplifying  assumption  of  the  previous 
section,  these  reference  frames  become  identical.  Only  F/  will  be  used  here.  Thus  let  F, 
be  an  Earth-fixed  inertial  reference  frame  whose  origin  is  at  the  launch  site,  whose  x-axis 
points  along  the  projection  of  the  nominal  launch  trajectory  onto  the  Earth’s  surface  and 
whose  z-axis  is  nominally  downwards  (see  Fig.  1).  The  y-axis  follows  from  the  right 
hand  rule. 

A  second,  inertial  Earth- fixed  reference  frame  that  is  useful  is  the  launcher 
reference  frame  FL.  The  origin  of  FL  is  located  at  the  launch  site  with  the  x-axis  pointing 
in  the  launch  direction  and  the  z-axis  being  nominally  downwards  (see  Fig.  1).  The 
y-axis  follows  from  the  right  hand  rule. 

A  third  Earth- fixed  reference  frame  which  is  occasionally  required  is  that  of  a 
reference  frame  FT  placed  at  some  distance  from  the  launch  site.  This  reference  frame 
may  be  used  to  compute  the  rocket  aspect  angle  presentations  (e.g.  from  the  training 
ship).  Since  the  location  and  orientation  of  this  reference  will  depend  on  the  particular 
application,  it  is  defined  only  generally  in  Figure  1 . 

Since  the  aerodynamic  forces  are  most  conveniently  expressed  with  respect  to  the 
vehicle,  a  body- fixed  reference  frame  FB  will  also  be  used.  The  origin  of  FB  is  located  at 
the  vehicle  centre- of- mass.  In  vehicles  that  are  axisymmetric,  the  x-axis  is  on  the  axis  of 
symmetry  and  points  forward  through  the  nose.  Otherwise  the  x-axis  points  in  the 
nominal  launch  direction.  The  z-axis  is  nominally  downward,  while  the  y-axis  follows 
from  the  right  hand  rule.  This  reference  frame  and  some  associated  aerodynamic  angles 
are  shown  in  Figure  2. 

For  cases  where  the  rocket  vehicle  mass  characteristics  are  axisymmetric,  the 
aerodynamic  forces  are  independent  of  the  vehicle’s  roll  attitude,  and  the  thrust  forces 
are  axisymmetric,  the  body- fixed  reference  frame  need  not  spin  with  the  vehicle.  Thus  a 
reference  frame  FB'  is  defined  which  is  identical  to  FB  except  that  it  does  not  rotate  with 
the  vehicle  about  the  axis  of  symmetry.  Initially  FB  and  FB'  will  coincide. 

A  third  body- fixed  reference  frame  that  is  useful  in  specifying  the  vehicle’s  mass, 
inertia  and  configuration  characteristics  is  a  reference  frame  F*  whose  origin  is  located 
on  a  nose  datum  plane  on  the  vehicle.  If  the  vehicle  is  axisymmetric,  then  the  origin  of  F„ 
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is  on  the  axis  of  symmetry,  and  its  x-axis  points  towards  the  rear  of  the  vhicle  on  the  axis 
of  symmetry.  Otherwise,  the  origin  of  F«  may  be  any  convenient  location  on  the  nose 
datum  plane,  and  the  x-axis  nominally  points  towards  the  rear  of  the  vehicle.  The  z-axis 
is  nominally  downward  and  the  y-axis  follows  from  the  right  hand  rule  (see  Figure  2). 

In  the  following,  use  is  made  of  a  number  of  notation  conventions.  In  particular 
if  Lba  denotes  a  rotation  matrix  relating  the  components  of  a  vector  V  expressed  in 
F^fV/t)  to  the  components  of  the  same  vector  in  Fs(vfl),  then 


vfl  =  Lba  V* 


(2.2,1) 


The  following  definitions,  geometric  relationships,  and  matrices,  will  be 
employed  in  the  presentation  of  the  equations  of  motion. 


1.  The  rotation  matrix  relating  F,  and  Fs: 

COS  0b  COS  4/fl 

cosOsSini^a 

-  sin0B 

Lb/  — 

sintasindBCOsi^B 
-  cos^sintps 

sin^aSinOasiny/fl 
+  COs4bCOS4>b 

sin<t>aCOS0a 

cos4Bsin0flcos4/fl 
+  sin^Bsiny/a 

cos4Bsin0Bsin4>B 
—  sin^aCOsifiB 

COS 4b cos 0B 

(2.2,2a) 


(2.2,2b) 


where  4b,  dB,  and  yB  are  the  Euler  angles  defined  by  Etkin  (Reference  3). 

The  rotation  matrix  relating  F,  and  F«'  follows  from  (2.2,2a)  by  substituting  yjB' , 
Qa’  and  4b'  for  ipB,  6g  and  4  s  respectively,  and  will  be  denoted  Lb'j. 


2.  The  rotation  matrix  relating  Ft 

and  F,: 

cos08o 

0 

-  sin0Bo 

L/,  = 

0 

1 

0 

(2.2,3a) 

sin0fl„ 

0 

COS0bo 

=  [1 

LlJ 

(2.2,3b) 
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3.  The  rotation  matrix  relating  Fr  and  F,: 

(2.2,4a) 

=  [lr/J  (2.2,4b) 

4.  The  angular  velocity  of  the  vehicle  with  respect  to  F/  written  as  components  in  F*  and 
F«' : 

oiB  =  (Pa,  Q«,  r*)r  (2.2,5a) 

a,*'  =  (0,  qfl',  Xs'Y  (2.2,5b) 


coster  sinter  0 

-  sim^r  coster  0 

0  0  1 


5.  The  angular  rate  cross-product  matrices  (Reference  3)  in  F*  and  F*’ : 


0 

q» 

to® 

= 

rB 

0 

-Pa 

(2.2,6a) 

~  Qs 

p« 

0 

0 

-  Ta' 

Qa' 

>\i 

O)®' 

= 

Ta' 

0 

0 

(2.2,6b) 

-q*' 

0 

0 

The  airspeed  vector  of  the 

:  vehicle  written 

as  components  in  Fa  and  Fa' 

: 

V* 

=  (Ua,  V«,  Wa)r 

(2.2,7a) 

ya* 

=  (Ua 

>»  V»',  Wa') 

r 

(2.2,7b) 

7.  The  groundspeed  vector  of  the  vehicle  with  respect  to  F,  written  as  components  in 
F*,  Fa'  and  F,: 


I 
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V®  =  (uflt.,  \be,  W Be)t 

(2.2,8a) 

V®'  =  (Ufl£,  vfl£\  wfl£.')r 

(2.2,8b) 

* 

Vr  =  (x„  y„  z,)T 

(2.2,8c) 

8.  The  aerodynamic  angles  (see  Figure  2): 


a 

=  arctan  (wfl/ua) 

(2.2,9) 

P 

=  arctan  (va/V„) 

(2.2,10) 

where 

=  Vu*1  +  Wfl1 

v„ 

(2.2,11) 

V 

=  V^u,,2  +  V„2  +  W  B2 

(2.2,12) 

Here  a  is  the  angle  of  attack  of  the  x-axis  of  F»,  /?  is  the  sideslip  angle  of  the  x-axis 
of  F„,  V  is  the  airspeed,  and  V„  is  the  magnitude  of  the  airspeed  vector  component  along 
the  x-z  plane  of  F». 

Analogous  angles  to  a  and  p  may  be  written  in  terms  of  Fa '  components  by  direct 
substitution  of  F„'  quantities  for  F„  quantities. 

9.  The  geometric  relationships  (see  Figure  2): 


ufl  =  Vcos^cosa 

(2.2,13a) 

va  =  Vsin  P 

(2.2,13b) 

ws  =  Vcos^sino 

(2.2,13c) 

The  wind  velocity  with  respect  to  F,  written  as  components  in  F„ 

F»,  and  Fb'  : 

£ 

M 

£ 

II 

(2.2,14a) 

w®  *  (Ub,,  v«g,  wBg)r 

(2.2,14b) 

w®'  =  (uv  Vb/,  wBg')r 

(2.2,14c) 
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11.  The  acceleration  due  to  gravity  written  as  components  in  F,: 


g'  =  (0,  0,  g)r  (2.2,15) 

12.  The  aerodynamic  forces  (not  including  thrust  forces)  written  as  components  in  F* 
and  F : 


A"  =  (XAb,Ya.,Z  Ab)t  (2.2,16a) 

A®'  =  (XAb,YAb',YAb')t  (2.2,16b) 

13.  The  aerodynamic  moments  (not  including  thrust  moments)  written  as  components 
in  F*  and  F*' : 


M®  =  (LAatMABt^AB)T  (2.2,17a) 

M®'  =  (LAb,  MAb',MAb')t  (2.2,17b) 

14.  The  inertia  matrix  of  the  vehicle  with  respect  to  its  centre -of- mass  expressed  in  F„ 
(see  Etkin,  Reference  3)  and  F*' : 


I® 

lXX 

-  Is 

Ixy 

-I® 

1*  = 

-I® 

l*y 

I® 

lyy 

-I® 

(2.2,18a) 

-i® 

*xz 

-I8 

i® 

lxx 

0 

0 

Is'  = 

0 

I®' 

lyy 

0 

(2.2,18b) 

0 

0 

I®' 

*yy _ 

15.  The  total  thrust  forces  written  as  components  in  F*  and  FB' : 


TB  =  (XrB,YrB,ZrBy  (2.2,19a) 

I®'  =  (Xr„  Yr.',  Ytb'V  (2.2,19b) 
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16.  The  total  thrust  moments  as  components  in  F„  and  F»' : 

=  (LTfl,  Mrg,  NrB)T  (2.2,20a) 

M V  =  (LTB,MTB',MrB’)T  (2.2,20b) 

2.3  Newton-Euler  Development  of  the  General  Equations  of  Motion 

Newton-Euler  techniques  begin  with  the  fundamental  equations  (Reference  3) 


F 

=  ma 

(2.3,1) 

=  M 

(2.3,2) 

a  is  the  acceleration  vector  of  the  body  centre -of- mass  relative  to  an  inertial  reference 
frame,  h  is  the  angular  momentum  of  the  body  about  its  centre-of-mass,  _F  is  the 
external  force  vector  acting  at  the  centre-of-mass  and  M  is  the  external  moment  vector 
about  the  centre-of-mass.  JF  may  be  written 

F  =  mV,  (2.3,3) 

where  V£  is  the  velocity  vector  of  the  vehicle  with  respect  to  F,  and  m  is  mass  of  the 
vehicle.  An  expression  for  follows  from  the  fundamental  relationship 

h  =  /  [  r  x  r  ]  dm  (2.3,4a) 

^  man  — *  * 

or 

h  =  /  [rxr  +  rx(co*xr)]dm  (2.3,4b) 

*  mass  *  *  *  * 

where  r  is  the  position  vector  of  an  element  of  mass  dm  of  the  body  with  respect  to  its 
centre-of-mass  (see  Figure  3),  co  *  is  the  angular  velocity  vector  of  F„  with  respect  to  F„ 
when  applied  to  a  vector  represents  rate  of  change  with  respect  to  F/  and  ‘o’  when 
applied  to  a  vector  represents  rate  of  change  with  respect  to  F*  (see  Reference  4  for  a 
more  thorough  discussion  of  vector  differentiation).  Equation  (2.3,4b)  may  be  written 
in  matrix  notation  as  (replacing  x  by  -  r*  x  cj*  and  dropping  the  subscript  ‘B’ 
on  co b  for  the  sake  of  brevity)* 

*  Superscripts  on  matrix  quantities  refer  to  the  reference  frame  in  which  the  components  of  the 
matrix  are  expressed.  Overscore  refers  to  the  matrix  equivalent  of  the  vector  cross-product. 


UNCLASSIFIED 


UNCLASSIFIED 


/9 


But 


ha 
r  8 


J  [i8r8-i8?8u>*]dm 


0 


(2.3.5) 

(2.3.6) 


for  a  rigid  body.  Since  o)“  is  a  constant  with  respect  to  the  integration  in  (2.3,5),  it 
follows  that 


h*  =  Xs  co 8  (2.3,7) 

where 

Is  =  -/  (2.3,8) 

moss 

X®  is,  by  convention,  given  by  (2.2,18a). 

The  externally  applied  force  £  is  made  up  of  an  aerodynamic  component  A,  a 
thrust  component  T  and  a  gravitational  component  mg  such  that 

F  =  A  +  T  +  mg  (2.3,9) 

Substituting  (2.3,9)  into  (2.3,1),  the  vector  force  equation  becomes 


mVf  =  A  +  T  +  mg  (2.3,10) 

The  externally  applied  moment  M  is  made  up  of  an  aerodynamic  component 
and  a  thrust  component  Mr  such  that 

M  =  M,  +  MT  (2.3;  1 1 ) 

Substituting  (2.3,11)  into  (2.3,2)  yields 

^  =  M„  +  Mr  (2.3,12) 

Other  than  the  gravitational  force,  the  dominant  forces  and  moments  acting  on 
the  aircraft  are  due  to  aerodynamic  causes  and  are  largely  determined  by  its  orientation 
and  configuration.  It  is  accordingly  advantageous  to  write  the  matrix  equations  of 
motion  with  respect  to  a  body- fixed  reference  frame.  This  reference  frame  is  chosen  to 
be  F*.  Furthermore,  this  choice  does  not  introduce  any  gravitational  moments  since  the 
origin  of  F«  and  the  centre-of-mass  of  the  vehicle  coincide. 


UNCLASSIFIED 


UNCLASSIFIED 


/IC 


Thus  the  matrix  force  and  moment  equations  become 

m(V|  +  a^Y!)  =  A*  +  T*  +  mUg' 
and 

hfl  +  hfl  =  M5  +  M? 

Substituting  for  h*  from  (2.3,7),  the  moment  equation  becomes 
\B6>*+  1  *  o)B  +  B  (oB  =  M5  +  M? 

Writing  out  the  equation  (2.3,13)  in  scalar  form  yields 

m(uSE  +  q^Wsg  -  rflVflE)  =  XAa  +  XTb  +  mgl«,13 

m(v„E  +  rBuflE  -  p*Wj,E)  =  Y**  +  YTa  +  mgl*/23 

m(wBE  +  p*v*  -  q0Ui>£)  =  ZAa  +  ZTa  +  mgl*^ 

for  the  force  equations,  and 

IU>b  -  I?„q*  -  I«rfl  +  If,p*  -  I*q*  -  i®r*  +  Iy,(ra  -  q l) 
+  (If,  -  Ify)rBq*  +  IfyrBpfl  -  If.qaP*  =  LAb  +  LTb 

-  Kf> B  +  I-q*  -  If,  fa  -  I®  Pa  +  ifyqa  -  if.r»  +  If,(Pa  -  r*) 

+  (I®  -  If.)p*rB  +  If, pB q B  -  I®  r*q*  =  M«B  +  Mr* 

-  I®  Pa  -  If,  qa  +  If,  fa  -  If,  Pa  -  If.qa  +  If.  Ta  +  I®(qB  -  pi) 

+  (I®  -  I®)Paqfl  +  If.qaTa  -  If.PaTa  =  N4*  +  Nr* 


(2.3.13) 

(2.3.14) 

(2.3.15) 

(2.3,16a) 

(2.3,16b) 

(2.3,16c) 

(3.2,17a) 

(3.2,17b) 

(3.2,17c) 


for  the  moment  equations. 

Kinematic  equations  are  also  required  for  the  linear  and  rotational  position  of  the 
aircraft.  The  linear  position  equations  follow  from 


Vi  =  L«  Vf 


(2.3,18) 
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The  rotational  position  equations  are  the  Euler  angle  rate  equations  and  are 
derived  in  Etkin  (Reference  3).  The  resulting  scalar  kinematic  equations  of  motion  are 
thus  seen  to  be 


X/  =  Ib/jjUbj.  +  lB/j,vS£  +  1«/31  wfl£  (3.2,19a) 

Y/  =  ls/I2U*£  +  1b/22V»£  +  l«/32W«f  (3.2,19b) 

Z/  =  ls/|3  Ujj£  +  \bi13Vbe  +  1«/33Wb£  (3.2,19c) 

for  linear  position,  and 

4>s  =  p B  +  q«sin<}>Btan0B  +  r«  cos<j>fl  tan0B  (2.3,20a) 

QB  =  qBcos<f>B  -  rflsin<t>B  (2.3,20b) 

Vb  =  [qssin<t»B  +  rBcos<t>fl]sec0B  (2.3,20c) 

for  angular  position. 


It  should  be  stressed  that  the  variables  (uBe,  vBe,  wBe)  in  equations  (2.3,16), 

(2.3,17)  and  (2.3,19)  are  the  body-axes  components  of  the  vehicle’s  ground  velocity 

vector.  This  is  not  the  same  as  the  equations  developed  in  Reference  2  where  airspeed 

vector  components  in  body- axes  are  used.  Also,  no  assumptions  have  been  made,  up  to  l 

this  point,  regarding  vehicle  planes  of  symmetry  and  the  symmetry  of  the  thrust  and 

aerodynamic  forces  and  moments.  Finally,  it  should  be  noted  that  no  assumptions  have 

been  made  about  the  mass  and  inertia  characteristics  of  the  vehicle,  i.e.  in  general 

I 

m  #  0 
and 

if,  *  0 

This,  too,  is  an  added  feature  not  present  in  the  equations  developed  in  I 

Reference  2. 


I 
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A  number  of  simplifications  may  be  added  to  these  equations  if  certain  symmetry 
conditions  are  satisfied.  If  the  vehicle  has  mass  symmetry  about  the  xy  and  xz  planes, 
then 

If,  =  If,  =  If,  =  0  (2.3,21) 

If  the  vehicle  mass  characteristics  are  also  axisymmetric,  then  in  addition  to 
(2.3,21)  we  also  have 

If,  =  If,  (2.3,22) 

Applying  assumption  (2.3,21)  to  the  moment  equations  (2.3,17a)  through  to 
(2.3,17c)  results  in  the  simplified  set  of  equations 

pa  =  [-if, Pa  -  (If,  -  If,)rflqa  +  L Ag  +  LrJ / Ifc  (2.3,23a) 

q»  =  l-  If,Qa  -  (If,  -  If,)Parfl  +  MAb  +  MrJ/I*  (2.3,23b) 

fa  =  [-  if,ra  -  (If,  -  If,)Paqa  +  NAg  +  NrjJ / If,  (2.3, 23c) 

Applying  the  axisymmetry  assumption  (2.3,22)  to  those  equations  simplifies 
(2.3,23a)  even  further  by  eliminating  the  (If,  -  If,)raq„  term,  i.e. 


Pa  =  [  —  Ifc  Pa  +  LAb  +  LrJ/I*  (2.3,24) 

If  we  now  make  the  assumptions  that  the  thrust  forces  are  axisymmetric,  that  the 
vehicle  mass  characteristics  are  axisymmetric1,  and  that  the  vehicle  aerodynamic 
characteristics  are  independent  of  the  vehicle’s  roll  orientation,  then  we  may  take 
advantage  of  the  simplifications  that  will  result  to  the  equations  of  motion  by  expressing 
them  in  the  reference  frame  Ffl'  rather  than  in  Fa.  Recall  that  in  Section  2.2  we  defined 
Fb'  as  being  identical  to  FB  except  that  it  does  not  rotate  with  the  vehicle  about  the  axis 
of  symmetry.  As  a  result,  if 

<vB'  =  psj^a  +  q«  +  rijcfl  (2.3,25) 

is  the  angular  velocity  vector  of  Ffl'  relative  to  the  inertial  reference  frame  F,  expressed 
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as  components  in  FB',  then  the  latter  condition  implies  that 

pB'  =  0  (2.3,26) 

for  all  t. 

Also  in  general  we  will  have 

qi  #  qB  (2.3,27a) 

r b  ^  r*  (2.3,27b) 

For  the  purpose  of  generating  the  aerodynamic  forces,  however,  qB  and  rB  may 
be  treated  as  being  interchangeable  with  qB  and  rB  respectively  because  of  the 
assumption  that  the  aerodynamic  forces  are  independent  of  the  roll  orientation. 
However,  the  resulting  aerodynamic  forces  are  now  expressed  in  FB  axes  rather  than  in 
Fb  axes. 

There  are  also  aerodynamic  forces  that  are  generated  due  to  the  vehicle’s  rolling 
angular  velocity  pB.  The  latter  is  not  identical  to  pB,  and  thus  an  equation  of  motion  for 
pB  will  still  have  to  be  retained. 

Formally  the  equations  of  motion  in  FB  may  be  obtained  from  the  equations  of 
motion  in  FB  by  setting  pB  =  0  in  all  equations  except  the  pB  equation,  replacing 
(u*£,  v„£,  wBj.)  with  (u„E,  \ZE,  wB£),  (qB,  rB)  with  (q;,  r^),  (<t»B,  M>b)  with 
(♦;,  vi),  (MAa,  NAb)  with (M;fl,  M;b),  (Mrfl,  NTb)  with (m;b,  m;8),  (X*b,  YAb>  ZAb) 
with  (X„B,  YAb,  Y\b),  (XrB,  Yt-.,  Z tb)  with  ( XT„ ,  YrB,  Y£B),  and  applying  the  assumptions 
discussed  previously.  Finally,  the  pB  equation  from  the  FB  equations  is  retained. 

The  resulting  equations  of  motion  are  as  follows: 
m(uBr  +  qiwBf  -  r;v;t)  =  XAb  +  XTb  +  mgU/l3  (2.3,28a) 

m(v;E  +  rBuBf)  =  Y;fl  +  Y;b  +  mgU/23  (2.3,28b) 

m(w,E  -  qBu„£)  =  Y;b  +  Y^b  +  mgU,33  (2.3,28c) 
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P» 

- 

(-  IS,  Pa  +  LAb  +  LrJ/I® 

(2.3,29a) 

q  b 

- 

[-  i?,qi  +  m;b  +  M^]/I* 

(2.3,29b) 

r k 

= 

[  -  I®  r«  +  M;,  +  Mtb]  /  Ify 

(2.3,29c) 

x, 

= 

lfl/,,  +  1  fl/2 ,  +  1  a/3 ,  Wfl£ 

(3.2,30a) 

y> 

= 

l«uUic  +  lfl/22  Vfl£  +  lfl/32  Wfl£ 

(3.2,30b) 

z  / 

= 

1b/,3  uBr  +  U/23  vb£  +  1«/33Wb£ 

(3.2,30c) 

k 

= 

qisin^stanSi  +  ricos^tanfl* 

(2.3,31a) 

% 

= 

q^cos^i  -  r*sin$i 

(2.3,31b) 

Vb 

= 

[qflSimj»a  +  ri  costal  sec0fl 

(2.3,31c) 

2.4  Aerodynamic  Model 

The  equations  of  motion  developed  in  the  previous  section  contain  terms 
(e.g.  X«)  that  represent  the  aerodynamic  forces  acting  on  the  vehicle.  In  this  section 
these  terms  are  defined  as  functions  of  the  vehicle’s  state. 

Although  more  sophisticated  techniques  are  available  (see  the  discussion  in 
Reference  2),  for  the  purposes  of  rigid  body  six  degree- of- freedom  simulation,  it  is 
usually  quite  adequate  to  use  a  quasisteady  aerodynamic  model  based  on  Bryan’s 
aerodynamic  derivative  technique. 

No  attempt  has  been  made  here  to  generalize  this  model  so  that  it  applies  equally 
well  to  all  types  of  air  vehicles.  Rather,  its  form  has  been  simplified  so  that  it  is  suitable 
for  use  only  with  free-flight,  ballistic,  rocket- boosted  vehicles. 

The  resulting  model  expressed  as  aerodynamic  force  and  moment  components  in 
F»  is  summarized  below.  No  attempt  is  made  to  rationalize  this  model  other  than  to  state 
that  its  use  has  resulted  in  predicted  trajectories  that  are  in  good  agreement  with 
measured  flight  characteristics  (see,  e.g.,  Reference  1)  of  CRV7/BATS  and  ROBOT-9 
vehicles. 
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The  aerodynamic  forces  are  specified  by 

\ab  =  ~  CdQdS  (2.4, 1  a) 

=  CypP  QoS  +  qD  S(Cya/1„  +  Cyoyjn)ps,udo  (2.4,1b) 

Z.4B  =  C,aoqDS  +  q  d  S  ( C ,  3/,,,  +  CIoyjn)p3tUdo  (2.4,1c) 

The  various  quantities  in  these  equations  are  defined  in  the  notation.  It  is 
important  to  note  that  qD  is  the  dynamic  pressure  given  by 

qD  =  ViqW1  (2.4,2) 

where  q  is  the  air  density  and  V  is  the  airspeed  given  by 

V  =  (u  *  +  v£  +  w  I)1'2  (2.4,3) 


Here  (u„,  vfl,  wfl)  are  the  airspeed  vector  components  expressed  in  F«,  and  in  the 
presence  of  nonzero  wind  conditions  will  not  be  t*-*  same  as  (uBjr,  v*E,  w„£).  Rather,  they 
will  be  related  to  the  wind  velocity  vector  components  in  F*  [(uBf,  \Bf,  wBt)]  through  the 
relationships 


V,  =  V  +  W 


(2.4,5) 


i.e.  the  ground  velocity  vector  V£  equals  the  airspeed  vector  V  plus  the  wind  velocity 
vector  W. 
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CD  is  the  nondimensional  drag  coefficient,  Cyp  is  the  aerodynamic  derivative 
relating  y-force  due  to  sideslip  angle  ft,  C,a  is  the  aerodynamic  derivative  relating  z- force 
due  to  angle  of  attack  a,  and  S  is  a  reference  area  that  is  usually  the  fuselage  cross- 
sectional  area  for  ballistic  vehicles.  CD,  Cyp,  and  C,a  may,  in  general,  be  Mach  number 
and  Reynolds  number  dependent,  although  here  they  are  considered  to  be  only  Mach 


number  dependent,  a  and  p  are  given  with  respect  to  the  x-axis 
geometric  considerations  may  be  shown  to  be  (see  Figure  2) 

of  Fb,  and  from 

a 

=  arctan(w„/uB) 

(2.4,6a) 

P 

where 

=  arctan(vfl/ V„) 

(2.4,6b) 

v„ 

=  (u2b  +  w2)1'2 

(2.4,7) 

Differential  equations  may  be  obtained  for  a  and  (i  by  differentiating  (2.4,6a)  and 
(2.4,6b)  with  respect  to  time,  with  the  results 

a  =  [w„/V  -  uB  wB/V2]coszcr  (2.4,8a) 

and 

P  =  [V„vB  -  vH(ufl uB  +  wBwB)V„]/ (V^cos/J)  (2.4,8b) 

The  last  terms  on  the  right  hand  sides  of  (2.4,1b)  and  (2.4,1c)  are  pseudo  fin 
terms  and  are  included  to  permit  modeling  of  aerodynamic  asymmetries  (e.g.  due  to 
production  tolerances).  They  do  not  include  any  of  the  effects  produced  by  the  vehicle’s 
nominal  fin  configuration.  The  latter  have  already  been  included  in  CD,  C yp  and  C,a. 

The  pseudo  fin  terms  are  defined  as  follows: 


r 

Cl^/in  ^ fin  Sin^/M 

(2.4,9a) 

c  = 

V'*o/(n 

—  Ciafin  6 fin  COS$//n 

(2.4,9b) 

Cj-a/it  = 

(2.4,10a) 

^ZQfin 

COS*|y,n 

(2.4,10b) 
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W/i.  =  wflcos<t»/1B  +  vBsin<|>/m  (2.4,11) 

afiH  =  arctan(w/iB/ufl)  (2.4,12) 

Here  CL(jfin  is  the  lift  slope  of  the  fin,  6fl„  is  the  cant  angle  of  the  fin,  is  the 
angular  cylindrical  coordinate  of  the  fin,  e/in  is  a  body- fin  interference  factor  and  a,in  is 
the  angle  of  attack  of  the  fin.  The  fin  geometry  coordinate  system  is  summarized  in 
Figure  4. 

The  aerodynamic  moment  expressions  are  defined  similarly  to  the  aerodynamic 
force  expressions,  as  follows: 


LAb  —  C,ppflqDSbV(2V)  +  C/^nd/iBqnSb 

Y ABnom  Zc*  ~  ZA Bnom  ^ cl 

=  Z4j,Bom(x«  -  X£<)  +  Cm#qBSqDb2/(2V) 

+  ZCJ  +  q  D  S  ( C ,  afin  Q fin  "t"  Ct0/ln)parudo{XaCfin  XCf) 

N4fl  =  -  Y^„0(n(xai  -  xc.)  +  CarrBSqDb2/ (2V) 

+  X-AgyCg  Qo  S  {Cyafin  +  o/i'n  )  psrudo  (  X  acjin  -  Xc,) 


(2.4,13a) 


(2.4,13b) 


(2.4,13c) 


In  these  expressions  Y4s(iom  and  ZABnom  are  given  by  (2.4,1b)  and  (2.4,1c)  without 
the  pseudo  fin  contributions,  (xct,  ycg,  zct)  are  the  coordinates  of  the  vehicle  centre-of- 
mass  in  the  vehicle  structural  reference  frame  F*  (see  Figure  2),  (x„c,  yac,  zoc)  are  the 
coordinates  of  the  vehicle  aerodynamic  centre  in  F„,  (xoc/in,  yac/(>i,  zoc/iK)  are  the 
coordinates  of  the  aerodynamic  centre  of  the  pseudo  fin  in  F*,  b  is  the  reference  length 
(usually  the  fuselage  diameter  for  ballistic  vehicles)  and  the  aerodynamic  derivatives  Clp, 
C,t iJiH,  Cmq,  Cia/in,  C„r,  Cya/in  are  defined  in  the  notation  list. 

The  aerodynamic  forces  and  moments  may  also  be  written  as  components  in  Fs  by 
making  an  identical  set  of  substitutions  into  (2.4,1)  and  (2.4,13)  as  used  in  converting  the 
equations  of  motion  written  in  Fs  to  those  written  in  FB'.  It  is  important  to  note  that 
certain  simplifications  result  because  of  the  underlying  assumptions  used  in  developing 
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the  F a'  equations,  i.e.  no  pseudo  fin  terms  may  be  included  and  y*  and  z’B  axes 
characteristics  are  identical. 

For  the  sake  of  brevity,  the  aerodynamic  force  and  moment  expressions  in  F*  will 
not  be  given  here. 

2.5  Mass  and  Moments  of  Inertia  Models 


The  equations  of  motion  have  been  written  so  that  variations  in  the  vehicle’s  mass 
and  inertia  characteristics  (due  to  rocket  motor  propellant  burn)  are  permitted. 
Component  methods  are  used  to  compute  the  total  vehicle  mass  and  moments  of  inertia. 
The  components  considered  are  the  vehicle  airframe,  the  vehicle  payload,  the  vehicle 
rocket  motors  less  propellant  and  the  rocket  motors’  propellant.  Of  these  components, 
only  the  propellant  characteristics  are  considered  to  be  variable  with  time.  Finally,  the 
assumption  has  been  made  that  the  payload,  the  rocket  motors,  and  the  propellant  are 
point  masses. 

Under  these  conditions  the  expressions  for  the  vehicle  mass  and  inertia 
characteristics  are  summarized  below.  These  equations  are  given  for  the  reference 
frame  Fs.  Position  coordinates  are  with  respect  to  the  vehicle  structural  reference 
frame  F*.  The  subscripts  used  to  reference  the  different  components  are  as  follows: 


1)  ‘em’  —  airframe  (empty) 

2)  ‘PL’  —  payload 

3)  ‘Me’  —  rocket  motors  less  propellant 

4)  ‘PR’  —  rocket  motor  propellant 

Nm  is  the  total  number  of  rocket  motors.  Other  variables  used  are  defined  precisely  in 
the  notation. 


The  expressions  for  the  mass  and  inertia  characteristics  are  given  by 
(AX|  =  x^  -  xc<,  Ay^  =  y^  -  yc„  and  so  forth) 


I 

m 

i 


Nm 

m<m  +  m„  +  Z  l(m„,),  +  (mm),] 


(2.5,1) 
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=  I®.*  +  m.*(Ay.J*  +  Az.2*)  +  mw(Ay^  +  AzPL) 

N„ 

+  .1  {(mOiKAy*.)?  +  (Az„.)2]  +  (m„),[(Ay„)?  +  (Az„)?]} 
=  I  ®>m  +  m.*(Ax?*  +  Az2*)  +  mPi(Axpi  +  AzPi) 

N„ 

+  . ? j  l(^*,)>l(AxM,)i  +  (Azi»,),]  +  (mPjt)<[( AxP#j)(  +  (AzPR),]} 

=  I®.m  +  m,*(Ax?*  +  Ay.2*)  +  mPL(Ax2PL  +  Ay’t) 

N„ 

+  .  ^  [(Axm,)<  +  (Ay*#,);]  +  (nipp),  [(  Axpr),  +  (AyP*),]} 


—  ~  m.„  Ax«*  Ay.*  mPi  AxPi  AyPI 
nm 

-  .  X  {(mM  ),( Axm,)/( Ay^.)/  +  (ni|.jl)i(Axi.ji)1(Ay/.u)/} 

i  =  1 


—  ~  m.iit  Ax.*  Az.*  —  m«  AxPi  AzPi 
N« 

~  {(mMe)<(^Xnf,),(AzM.),  -I-  ( mPJ, ) * ( AxPJt) , ( Azp„ ) , } 


—  ~  m.*  Ay.*  Az.*  —  vc\pl  Ay  pi  AzPl 
N« 

-  .X  {(mwJ.IAyMj^AzM.),  +  (mwrMAyroMAzp*),} 


{m.*x,*  +  mpix,.!  +  .X  [(m*.),(x*.),  +  (mPj,),(xPJ,)J}/m 


nm 

{m.*y.*  +  mPLyPL  +  .X,  [(nu.My*,),  +  (mpR),(ypR),]}/m 

i  =  l 


{m.*y,*  +  mptyPL  +  X  [(rTu.Mz*,^  +  (mpRMzpRLlJ/m 


UNCLASSIFIED 


UNCLASSIFIED 


/20 


Because  it  has  been  assumed  that  the  only  mass  changes  are  due  to  propellant 
burn,  in  these  expressions  the  only  time  variable  quantities  will  be  (mPR),,  (xPRi,  yPR(, 
z,,.).  If  we  further  assume  that  the  propellant  burns  in  such  a  way  that  the  centre- of- 
mass  of  the  propellant  of  a  given  rocket  motor  does  not  change  significantly  (e.g.  as 
would  be  the  case  in  rocket  motors  that  are  not  end  burners),  then  (xW()  yPRf,  zPR.)  are 
not  time  variable  and  only  (mPR),  need  be  considered.  The  latter  is  related  to  the  specific 
impulse  of  the  rocket  motor  through  the  relationship  (Reference  5) 

(mPR),  =  T,(t)/(liPl.g)  (2.5,4) 

where  g  is  the  acceleration  due  to  gravity,  T.(t)  is  the  thrust  of  the  i-th  rocket  motor  as  a 
function  of  time  t,  and  I,Pi  is  the  specific  impulse  of  the  i-th  motor. 

Under  these  assumptions  and  with  (2.5,4),  m,  the  moment  of  inertia  time 
derivatives,  and  (xcg,  ycg,  zcg)  may  be  readily  computed.  For  the  sake  of  brevity,  an 
exhaustive  set  of  equations  will  not  be  given.  Typically  we  have 

m  =  .1*  r-T,(t)/(I,p<g)]  (2.5,5) 

i  =  i 

It  =  -  2m„,(Ay„„ycg  +  Az,mzcg)  -2mP£(Ayycg  +  AzPizcg) 

N«  ... 

+  .  {  -  2(mM,).  [(AyM,).yc,  +  (AzM>),zc,]  +  (mPR),  (2.5,6a) 

[(AyPP)j  +  (Azpr),  ]  —  2  (mPR)i  [( AyPji)(  ycg  +  (AzPR  )(Zcg]} 

i?y  =  rn,m(xcgAy.m  +  Ax.myeg)  +  mPi(xctAyPt  +  yc*AxPI) 

N„ 

-  .  {  -  (mM€),[xct(AyM.)<  +  yct(AxM,),]  +  (mPR),  (2.5,6b) 

(  Axpr),  ( AyPR),  —  ( mPR)j  [xcg  ( AyPR),  +  ycg(AxPR),]} 

Nm 

xcg  =  {.Zj  I-T^O/d^JgKx™),  -  Xcgrh}m-  (2.5,7) 
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The  mass  and  moment  of  inertia  characteristic  formulations  expressed  in  the 
nonrotating  reference  frame  Ffl'  follow  from  the  Ffl  expressions  by  incorporating  the 
simplifying  assumptions  used  for  writing  the  Fs'  equations  of  motion  (see  Section  2.3). 
In  particular,  we  have 

If,  =  If,  =  If,  =  0  (2.5,8) 

and 

If,  =  If,  (2.5,9) 

Equations  (2.5,8)  and  (2.5,9)  are  just  the  result  of  the  mass  and  inertia 
axisymmetry  assumption  used  in  developing  the  Ffl'  equations  of  motion. 

For  the  sake  of  brevity,  the  expressions  for  the  F«'  mass  and  inertia  characteristics 
will  not  be  given  explicitly. 

2.6  Thrust  Characteristics 

In  Section  2.3  the  equations  of  motion  were  written  in  the  reference  F„  with  the 
thrust  forces  and  moments  written  generally  as  (Xrjg,  Yrs,  ZTg)  and  (LTb,  Mrs,  Nrfl) 
respectively.  In  this  section  these  terms  are  examined  in  more  detail. 

The  force  terms  (Xrfl,  YTg,  ZTg)  depend  on  the  time  domain  thrust  characteristics 
and  the  physical  location  and  orientation  of  the  rocket  motors.  This  data  must  be  known 
a  priori  to  the  simulation  and  is  provided  as  input  data  to  the  computer  program  in  the 
form  of  the  thrust  versus  time  look-up  tables.  The  transformations  used  are  summarized 
in  Appendix  2. 

The  thrust  moments  require  a  somewhat  more  detailed  examination.  They  are 
considered  to  consist  of  two  components: 

1)  A  moment  due  to  the  location  and  orientation  of  the  thrust  vector 
relative  to  F»  (see  Figure  4), 

2)  A  moment  induced  due  to  fixed  vanes  or  nozzle  grooves  onto  which 
the  exhaust  jet  impinges. 

Thus  we  have 

Lrs  =  (Lr,)«  +  (Lrfl)n,  (2.6,1a) 
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Mrfl  =  (Mr,),  +  (M,,)..  (2.6,1b) 

N  rfl  =  (NTb)c,  +  (Nr,)«  (2.6,1c) 

These  characteristics  are  rocket  motor  specific.  It  is  assumed  that  such  data  is 
available  for  the  rocket  motors  used  in  the  simulation.  It  then  follows  that  once  the 
orientation  and  location  of  the  rocket  motor  thrust  vectors  relative  to  the  vehicle  are 
specified,  enough  information  is  available  to  determine  (Lr,,  Mrj},  Nr„)  as  given  by 
(2.6,1a)  to  (2.6,1c)  (a  detailed  treatment  is  given  in  Appendix  2). 

An  assumption  that  has  tacitly  been  made  in  this  description  of  the  thrust  effects 
is  that  Coriolis  forces  and  moments  on  the  vehicle  generated  by  the  rocket  motor  exhaust 
are  negligible.  This  need  not  always  be  the  case,  particularly  for  the  moments,  if  the 
exhaust  mass  flow  rate  rh  and  the  exhaust  velocity  vector  relative  to  the  vehicle  V  *  are 
large.  However,  for  vehicles  in  the  class  of  CRV7/BATS  and  ROBOT-9  using  short 
burn  duration  70  mm  (2.75  inch)  rocket  motors,  these  effects  are  negligible  and  will  not 
be  considered  further  in  this  report. 

2.7  Vehicle  Kinematic  Restrictions  While  on  Launcher 

The  presence  of  the  launcher  during  the  initial  portion  of  the  flight  places  a 
number  of  kinematic  constraints  on  the  vehicle’s  motion.  This  section  considers  these 
constraints  for  a  rail  launcher  such  as  was  used  for  CRV7/BATS  and  ROBOT-9  (see 
Reference  1). 

The  basic  geometrical  quantities  are  defined  in  Figure  5.  The  equations  of  motion 
while  the  vehicle  is  on  the  rail  are  presented  for  the  following  assumptions: 

1)  The  vehicle  is  mechanically  constrained  from  tipping  backwards  or 
forwards  by  the  guide  T-bolt  until  the  T-bolt  clears  the  launch  rail 
(i.e.  the  vehicle  is  initially  constrained  to  move  along  the  x-axis  of 
reference  frame  FL).  In  the  case  of  CRV7/BATS  and  ROBOT-9 
there  is  also  the  launcher  cage  constraining  the  vehicle  for  part  of  its 
travel  on  the  launch  rail  (see  Reference  1). 

2)  The  quantity  sc  represents  the  distance  the  vehicle  must  move  in  the 
x-direction  of  F«  before  the  guide  T-bolt  clears  the  launch  rail.  The 


UNCLASSIFIED 


UNCLASSIFIED 


/23 


bolt  is  assumed  to  be  back  far  enough  on  the  vehicle  so  that  no 
significant  tip-off  may  occur  after  it  is  clear  of  the  rail  and  prior  to  the 
whole  vehicle  coming  clear. 

3)  The  vehicle  may  not  move  backwards  on  the  launcher  rail  (i.e.  \iB  is 
never  less  than  zero). 

Under  these  assumptions  it  follows  that  if  the  distance  that  the  vehicle  centre-of- 
mass  has  travelled  (s)  is  less  than  or  equal  to  sc,  then  the  vehicle  is  physically  constrained 
to  move  only  in  the  launch  rail  direction,  i.e.  for  s  <  sc  we  have 


Pa 

= 

O 

II 

II 

•> 

II 

«Q 
■  l— 

II 

ai 

•O” 

(2.7,1a) 

U« 

0 

(2.7,1b) 

VB 

= 

Va(0) 

(2.7,2a) 

w  B 

= 

w„(0) 

(2.7,2b) 

Pa 

= 

Pa(0) 

(2.7,2c) 

Qa 

= 

qa(0) 

(2.7, 2d) 

ra 

— 

ra(0) 

(2.1,2c) 

The  nonzero  conditions  (2.7,2a)  to  (2.7,2e)  allow  for  a  nonstationary  launcher, 
i.e.  as  would  be  the  case  for  a  launch  from  a  ship  in  linear  and  angular  motion. 

The  quantity  s  is  defined  precisely  as 

s  £  V(x,  -  x,.)1  +  (y,  -  y,o)J  +  z }  (2.7,3) 

where  (x„  y„  z,)  are  the  vehicle  centre-of-mass  coordinates  in  F,  and  (x,0,  y,o,  0)  are  the 
centre- of- mass  coordinates  when  the  vehicle  is  at  rest  on  the  launcher  prior  to  first  stage 
ignition. 

For  s  >  sG,  the  governing  equations  are  the  unconstrained  equations  of  motion 
developed  in  Section  2.3. 
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2.8  Wind  Model 

The  aerodynamic  model  presented  in  Section  2.4  includes  the  wind  velocities  in  F„ 
(F* ' )  as  (u„t,  vv  w„f)  [(ufl<,  w*f)J,  and  has  tacitly  assumed  that  there  is  no  variation 
of  the  wind  velocity  from  one  point  to  another.  This  is  equivalent  to  assuming  that  the 
wind  induced  aerodynamic  loads  are  determined  by  its  velocity  rate  of  change  acting  at 
the  centre- of- mass  of  the  vehicle,  an  assumption  referred  to  as  the  uniform-gust 
approximation  (References  2  and  3).  This  approximation  is  equivalent  to  assuming  that 
the  wind  velocity  spectral  content  significantly  affecting  the  vehicle  response  is  at 
wavelengths  that  are  greater  than  the  significant  vehicle  dimensions  (Reference  3).  This 
assumption  is  reasonable  when  considering  the  rigid  body  dynamic  response  of  flight 
vehicles,  particularly  for  smaller  vehicles  such  as  ROBOT-9  and  CRV7/BATS. 

The  wind  velocity  vector  components  relative  to  F,  are  most  conveniently 
expressed  as  components  in  F,,  i.e.  (Wlt  W2,  W3).  These  components  may  then  be 
related  to  the  wind  velocity  components  in  Ffl  with  the  rotation  matrix  Ls/  as  given  by 
(2.2,2a),  i.e. 

(«..,  vv  wfl<)r  =  L*,(Wt,  W2,  W 3)T  (2.8,1) 

The  simulation  package  has  provisions  for  inputting  (W„  W2,  W3)  as  functions  of 
altitude.  This  allows  modeling  of  wind  velocity  atmospheric  boundary  layer  effects, 
vehicle  encounters  with  jetstream  regions,  and  so  forth.  As  well,  since  meteorological 
winds  aloft  data  is  usually  given  as  a  function  of  altitude,  simulation  of  measured  wind 
conditions  is  facilitated. 

2.9  Atmospheric  Conditions 

Since  the  ROBOT-9  and  CRV7/BATS  vehicles  have  the  capability  to  achieve 
altitudes  well  above  9000  m  (30,000  ft),  an  atmospheric  model  is  required  that  takes  into 
account  variations  in  density  (g),  temperature  (T^),  pressure  (p>»),  and  the  speed  of 
sound  (a)  as  a  function  of  altitude  above  sea  level  (h4ii). 

The  models  used  are  based  on  the  U.S.  standard  atmosphere  (1962),  as  is  common 
practice  in  aeronautical  engineering,  and  are  valid  within  the  troposphere,  i.e.  for 
lust  <  11,100  m  (36,000  ft)  (see  Reference  6).  They  are  given  by 


UNCLASSIFIED 


UNCLASSIFIED 


/25 


Ta  =  288.15  -  0.0065  R  (2.9,1) 

PA  =  101300.(T/4/288.15)S2SS  (2.9,2) 

a  =  20.0463  (2.9,3) 

6  =  0.00348454  p^/T^  (2.9,4) 

where 

R  =  xE/  (h^5i  +  r£)  (2.9,5) 


Ty4  is  in  degrees  Kelvin,  p*  is  in  Pascals,  a  is  in  meters  per  second,  q  is  in  kilograms  per 
meter  cubed,  hASL  is  in  meters,  and  r£  is  the  Earth’s  radius  to  the  sea  level  datum, 

r.  =  6.3567658  x  106  m  (2.0855531  x  107  ft). 

From  the  factor  R  it  is  also  convenient  to  compute  the  variation  of  the 
acceleration  due  to  gravity  as  a  function  of  altitude,  i.e. 

g  =  goR2  (2.9,6) 

where  gn  =  9.80667  m/s2  (32.1741  f/s2). 

Provision  has  been  made  in  the  simulation  package  to  vary  the  temperature  and 
pressure  (and  thus  the  density)  from  the  standard  values  by  allowing  altitude  dependent 
per  cent  deviations  from  standard  conditions. 

2.10  Aspect  Angle  Equations 

For  target  and  flight  test  applications,  it  is  frequently  necessary  that  the  vehicle’s 
aspect  azimuth  (4,*)  and  elevation  (|£)  angles  be  known  with  respect  to  an  observer  at  Fr 
(see  Figure  l).  This  section  presents  equations  for  $A  and  4e  in  terms  of  the  location  and 
orientation  of  Fr  relative  to  that  of  F,. 

It  is  assumed  that  the  x-y  planes  of  F,  and  Fr  are  parallel,  i.e.  that  F,  may  be 
rotated  to  Fr  through  a  rotation  \^T  about  the  z-axis  of  F,. 
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Let  the  vector  position  of  Fr  relative  to  F/  be  Rr/,  and  that  of  the  vehicle  centre  - 
of-mass  relative  to  the  origin  of  F/  be  R  (see  Figure  6).  It  follows  that  the  vector 
position  of  the  vehicle  relative  to  Fr  is  given  by 

RT  =  R  -  R„  (2.10,1) 

or  in  matrix  notation 

R?  =  L  77  R 7  -  L  77  R  77  (2.10,2) 

where  L77  is  the  rotation  matrix  rotating  vector  components  in  F/  to  components  in  Fr, 


and  is  given  by  (2.2,4a). 

Equations  (2.10,2)  may  be  written  in  scalar  form  as 
Xr  =  (x,  -  X77)  cosier  +  (y /  -  yr/)sini^r  (2.10,3a) 

yT  =  -  (x,  -  xr/)sini^r  +  (y/  -  yr/) cosier  (2.10,3b) 

zT  -  z,  -  Z77  (2.10,3c) 

From  the  definition  of  and  in  Figure  1,  it  follows  that 

=  arctan  (yT/xr)  (2.10,4a) 

|£  =  arctan  (- zr/xr)  (2.10,4b) 


3.  BALSIM  SOFTWARE  DESCRIPTION  —  GENERAL 

The  dynamic  model  described  in  the  previous  chapter  has  been  implemented  in  the 
BALSIM  simulation  package.  All  coding  was  carried  out  using  IBM  FORTRAN  for  the 
H-extended  compiler.  The  package  has  been  debugged  and  tested  on  the  IBM  3033 
computer,  and  has  been  used  to  predict  the  dynamic  characteristics  of  the  CRV7/BATS 
and  ROBOT-9  vehicles. 


UNCLASSIFIED 


UNCLASSIFIED 


/27 


The  software  is  currently  being  installed  on  VAX1 1/780  and  Honeywell 
DPS-8/70C  computers. 

The  software  consists  of  a  MAIN  program  plus  nine  subroutines  making  up 
approximately  885  FORTRAN  source  statements.  There  are  no  subroutines  or 
functions,  other  than  these,  that  are  not  available  in  standard  FORTRAN  on-line 
libraries. 

The  software  userbook  is  given  in  Appendix  1  with  a  source  language  listing  of  the 
package. 

3.1  Software  Capabilities 

The  dynamic  model  implemented  with  the  BALSIM  package  has  already  been 
discussed  in  detail  in  the  previous  chapter.  Its  limitations  will  not  be  considered  further 
here. 


The  package  was  developed  with  the  objective  of  providing  a  convenient  basis  for 
inputting  the  characteristics  of  multistaged  rocket  vehicles  and  predicting  their  dynamic 
rigid  body  characteristics.  By  appropriately  modifying  the  input  data  set,  it  provides  for 

1)  nominal  and  off-nominal  vehicle  mass,  inertia  and  thrust 
characteristics, 

2)  different  motor  types, 

3)  Mach  number  dependent  aerodynamic  characteristics, 

4)  structural  production  tolerances, 

5)  system  failures  (e.g.  stage  and  fin  failures), 

6)  moving  launchers, 

7)  user  specified  initial  conditions, 

8)  user  specified  payload  characteristics, 

9)  tabular  output  in  either  metric  or  English  units,  and 

10)  multiple  case  runs. 

As  well,  with  minor  software  modification,  response  calculations  may  be  stored 
on  disk  for  subsequent  use  with  other  software  (e.g.  plotting  software). 
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3.2  Software  Limitations 

In  its  current  form  the  software  is  not  intended  for  use  in  the  following  types  of 
simulations: 

1)  Ballistic  rocket  vehicles  with  control  surfaces. 

2)  Winged  flight  vehicles. 

3)  Nonrigid  vehicles. 

4)  Simulations  where  the  inertial,  flat  Earth  approximations  are  invalid. 

5)  Vehicles  where  the  staging  process  involves  physically  releasing  rocket 
motor  stages. 

Limitations  (1)  and  (5)  may  be  removed  with  relatively  minor  alterations  to  the 
dynamic  model  of  Chapter  2  with  corresponding  changes  to  the  software. 

3.3  Numerical  Integration  Algorithm 

The  numerical  integration  algorithm  used  to  solve  the  system  of  ordinary 
differential  equations  describing  the  vehicle’s  dynamics  is  a  fixed  step-size,  fourth  order 
Runge-Kutta  method  (see  Reference  7).  Provision  has  been  made  for  the  specification 
of  two  step  sizes,  one  for  use  during  rocket  motor  burns,  and  the  other  for  use  during 
coasting  flight.  The  latter  technique  was  found  to  considerably  reduce  CPU  time  in 
certain  simulations. 

3.4  Software  Testing  and  Execution  Times 

The  BALSIM  package  has  been  used  extensively  to  predict  the  performance  and 
dynamic  characteristics  of  the  CRV7/BATS  and  ROBOT-9  vehicles.  These  predictions 
have  been  used  to  define  the  nominal,  dispersion  and  safety-envelope  characteristics  (see 
References  8  and  9)  of  these  vehicles. 

Flight  test  data  obtained  early  in  the  development  of  CRV7/BATS  (see 
Reference  1)  was  used  to  update  and  validate  the  aerodynamic  model  that  has  been 
employed.  More  recent  comparisons  with  flight  test  data  have  also  proven  to  be 
satisfactory  (also  see  Reference  1). 
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BALSIM  predictions  have  also  been  evaluated  for  consistency  by  comparing 
results  obtained  using  the  equations  of  motion  written  in  F«  with  those  written  in  F*' , 
with  satisfactory  results. 

The  execution  CPU  time  of  the  package  will  depend  on  the  computer  used,  on  the 
step  sizes  chosen,  and  on  the  duration  of  the  flight  time  simulated.  For  the  IBM  3033 
computer,  the  following  CPU  execution  times  were  observed  for  simulations  of  the 
ROBOT-9  vehicle  using  a  0.05  second  integration  step  size,  and  a  1.0  second  tabulated 
output  increment: 

1)  For  8  cases  averaging  108  simulated  flight  seconds  per  case,  the 
execution  CPU  time  was  343.4  seconds,  yielding  a  0.4  seconds  CPU 
execution  time  per  simulated  flight  second  ratio. 

2)  The  compilation  CPU  time  with  the  H-extended  compiler  was 
25  seconds. 

3)  The  linkage  editor  CPU  time  was  1.6  seconds. 

This  completes  the  general  description  of  the  BALSIM  software  package. 
Detailed  user  related  data  is  given  in  Appendix  1 . 

4.  SUMMARY 

Six  degree-of-freedom,  rigid  body  equations  of  motion  suitable  for  modeling  the 
dynamic  characteristics  of  multistaged,  free- flight,  ballistic  rockets  have  been  rigorously 
developed,  and  have  been  implemented  in  a  FORTRAN  software  package  called 
BALSIM.  This  package  allows  for  modeling  of  vehicle  thrust  and  structural 
asymmetries,  time-varying  mass  and  inertia  characteristics,  variable  wind  conditions, 
nonstandard  atmospheric  conditions,  stage  failures,  and  different  rocket  motor  types. 

The  BALSIM  package  has  been  successfully  used  to  predict  the  performance  and 
dynamic  characteristics  of  the  CRV7/BATS  and  ROBOT-9  vehicles  both  with  and 
without  moving  launchers.  It  will  be  adapted  for  use  with  the  VAX1 1/780  and 
Honeywell  DPS-8  computers  in  the  near  future. 
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DEFINITION  OF  F„,  FB',  FR,  o  AND  p 
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FIGURE  3 

ANGULAR  MOMENTUM  CONTRIBUTION  OF 
THE  MASS  ELEMENT  dm 
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FIGURE  4 

VEHICLE  CYLINDRICAL  COORDINATE  SYSTEM 
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s  =  distance  vehicle  has  moved 
along  launcher  rail  relative 
to  rest  position 


FIGURE  5 

VEHICLE  KINEMATIC  CONSTRAINTS  WHILE  ON 
LAUNCHER  GEOMETRY 
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ASPECT  ANGLE  GEOMETRY 
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BALSIM  USERBOOK 

This  appendix  contains  detailed  user  information  for  the  BALSIM  package.  A 
source  listing  of  the  IBM  H- Extended  FORTRAN  version  of  the  package  is  given  at  the 
end  of  the  Appendix. 

As  far  as  is  possible,  the  user  instructions  are  not  system  specific,  and  will  in 
general  apply  to  BALSIM’s  implementation  on  the  VAX1 1/780  and  the  Honeywell 
DPS-8/70C  computers  as  well  as  for  the  IBM  3033. 

General  background  information  on  BALSIM  is  contained  in  Chapter  3  of  the 
main  text  of  the  report.  The  equations  of  motion  that  are  used  are  developed  in  detail  in 
Chapter  2. 

Al.l  PROGRAM  UNITS 

The  BALSIM  package  consists  of  a  MAIN  program  plus  nine  subroutines.  These 
subroutines  are  as  follows: 

1)  INPUT  —  input  data  from  a  card  image  file  associated  with  unit 
number  NIN. 

2)  INPUS  —  modifies  input  data  as  required  for  the  NFIN  cases. 

3)  RUNK1 — performs  one  step  of  the  fourth  order  Runge-Kutta 
integration. 

4)  ATMOS  —  defines  U.S.  Standard  Atmosphere  (1962)  and 
acceleration  due  to  gravity  characteristics  as  a  function  of  altitude 
above  sea  level  (ALT). 

5)  WIND  —  defines  components  of  the  wind  velocity  vector  in  body- 
axes  as  a  function  of  altitude  above  sea  level  (ALT). 

6)  INTEQ  —  performs  a  linear  interpolation  of  arrays  TY1,  TY2,  TY3, 
whose  common,  possibly  unevenly  spaced  abscissae  are  stored  in  the 
array  TX. 
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7)  INTPL  —  performs  a  linear  interpolation  of  the  array  TY  whose 
abscissae  are  stored  in  the  array  TX. 

8)  OUTPU  —  carries  out  tabular  output  as  required  by  the  parameter 
NCOUT. 

9)  OUTAB  —  carries  out  summary  output  for  all  cases  run. 

A1.2  INPUT  DATA 

The  input  data  is  stored  on  a  card  image  file  associated  with  unit  number  NIN. 

The  input  data  organization  is  summarized  below,  and  may  also  be  obtained  from 
subroutine  INPUT. 


DATA  ITEM  FORMAT 

1.  NOROL,  NCOUT,  NIOUT  712 

2.  DELI,  DEL2,  TSTOP,  ALTST,  DELPT  8F10.2 

3.  SGUID,  SFRLN  8F10.2 

4.  XIO,  YIO,  ALTO  8F10.2 

5.  UO,  VO,  WO  8F10.2 

6.  THETO,  PSIO,  PHIO  8F10.2 

7.  PO,  QO,  RO  8F10.2 

8.  CLP,  CMQ,  CLDA,  DELTF,  B,  SREF  8F10.2 

9.  CLFIN,  DDFIN,  CPFIN,  PHFIN,  FIFAC  8F10.2 

10.  DXCP,  FXCN,  FXCD  8F10.2 

11.  NMACH  712 

12.  Repeat  NMACH  times: 

AMACH(I),  CDOT(I),  CNAT(I),  XCPT(I)  8F10.2 

13.  WTAPE,  XCGTP,  RCGTP,  PCGTP,  X1XTP,  XIYTP,  D1YZ  8F10.2 

14.  WPLD,  XCGPL,  RCGPL,  PCGPL  8F10.2 

15.  MTYPE  712 
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DATA  ITEM  (cont’d) 

FORMAT 

16. 

Repeat  MTYPE  times: 

WTE(l),  WTPR(l),  SPIMP(I),  CTHR(l),  CTHRT(I) 

8F10.2 

NTTHR(l) 

712 

TTHRT(1,1),  TTHR(l.l) 

8F10.2 

• 

• 

• 

• 

• 

• 

TTHRT(NTTHR(I),I),  TTHR(NTTHR(I),I) 

8F10.2 

17. 

NKICK 

712 

18. 

Repeat  NKICK  times: 

NTYPE(I),  NMOT(I) 

712 

TIG(I),  XTHR(I),  RTHR(I),  PTHR(I),  DTHR(I),  ATHR(I) 

8F10.2 

XEMPT(I),  XPROP(I),  REMPT(I),  RPROP(I),  PEMPT(I), 

PPROP(I) 

8F10.2 

TORQL(I),  TCCD(I),  TCCL(I),  TCCP(I) 

8F10.2 

19. 

NALTW 

712 

20. 

If  NALTW  <  0,  omit.  Otherwise  repeat  NALTW  times: 
ALTW(I),  WIXA(I),  WIYA(I),  WIZA(I) 

8F10.2 

21. 

NALTA 

22. 

If  NALTA  <  0,  omit.  Otherwise  repeat  NALTA  times: 
ALTA(I),  TMPR(I),  PRES(I) 

8F10.2 

23. 

NFIN 

712 

This  input  data  is  defined  as  follows: 

1.  NOROL  =  I 

a)  I  =  1  -  Ffl'  equations  of  motion  (axisymmetrical  mass, 
inertia,  aerodynamic  and  thrust  characteristics). 

b)  I  =  0  -  Fb  equations  of  motion  (axisymmetrical  mass, 
inertia  and  aerodynamic  characteristics). 

c)  I  =  -1  -  Fb  equations  of  motion  (general  case). 
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2.  NCOUT  =  IJ:  Tabular  output  control  parameter. 

a)  I  =  0,  J  >  0  —  Output  tables  1  to  J  +  1  inclusive  in 
imperial  units. 

b)  I  =  1 ,  J  >  0  —  Output  tables  1  to  J  +  1  inclusive  in 
MK.S  units. 

c)  I  <  0  —  Output  no  tables. 

3.  NIOUT  =  I:  Input  data  output  control  parameter. 

a)  1  =  0  —  No  listing  of  input  variables. 

b)  1  =  1  —  Listing  of  input  variables. 

4.  DELI:  Full  Runge-Kutta  integration  step  size  with  rocket  motors 
off  (seconds). 

5.  DEL2:  Full  Runge-Kutta  integration  step  size  with  rocket  motors 
on  (seconds). 

6.  TSTOP:  Maximum  permissable  simulation  time  (seconds). 

7.  ALTST:  Altitude  above  launch  point  at  which  the  simulation  will 
end  (ft). 

8.  DELPT:  Time  step  for  output  tables  (seconds),  tm„/DELPT  + 

1  «  200. 

9.  SGUID:  sc  (see  Figure  5,  ft). 

10.  SFRLN:  s„  (see  Figure  5,  ft). 

1 1 .  XIO,  YIO,  ALTO:  x,  (0),  y ,  (0),  hASL  (0)  (ft) 

a)  lw  (0)  is  the  altitude  above  sea  level  of  the  centre- of- 
mass  of  the  vehicle  while  at  rest  on  the  launcher. 

12.  UO,  VO,  WO:  u,,  (0),  vflf  (0),  w*r  (0)  (fps). 

13.  THETO,  PSIO,  PHIO:  Qa  (0),  tpB  (0),  (0)  (deg). 
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14.  PO,  VO,  RO:  p*  (0),  v,  (0),  .»  (0)  (deg/ second). 

15.  CLP:  C,p  (relative  to  the  diameter  and  cross-sectional  area  of  the 
fuselage,  1/deg). 

16.  CMQ:  Cm?  (relative  to  the  diameter  and  cross-sectional  area  of  the 
fuselage,  1/deg). 

17.  CLDA:  C, . ,  (relative  to  the  iiameter  and  cross-sectional  area  of  the 
fuselage,  1/deg). 

18.  DELTF:  dfim  (deg). 

19.  B:  Fuselage  diameter,  b  (ft). 

20.  SREF:  Fuselage  cross-sectional  area,  S  (ft2). 

21.  CLFIN:  CLg  pseudo  fin,  relative  to  S  (1/deg). 

22.  DDFIN:  (d/4„)  pseudo  (deg). 

23.  CPFIN:  (ft). 

24.  PHFIN :  angular  cylindrical  coordinate  of  aerodynamic  centre  of 

the  pseudo  fin  (see  Figure  4,  deg). 

25.  FIFAC:  e/in,  pseudo  fin  efficiency  factor,  >0  positive  fin,  <0 
negative  fin. 

26.  DXCP:  Vehicle  aerodynamic  centre  correction  term,  XCPT  = 

XCPT  +  DXCP  (positive  values  move  the  aerodynamic  centre 
further  back  on  the  vehicle,  ft). 

27.  FXCN:  Factor  adjusting  C*ra  data  (normally  FXCN  =  1),  i.e. 

Cit  =  C„0  *  FXCN. 

28.  FXCD:  Factor  adjusting  CD  data  (normally  FXCD  =  1),  i.e. 

d  =  Co  *  FXCD. 

29.  NMACH:  The  number  of  different  Mach  numbers  at  which 

aerodynamic  data  is  given  (NMACH  <  20). 
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30.  AMACH :  Mach  number  array  of  dimension  20. 

31.  CDOT:  Cd(M0)  array  of  dimension  NMACH. 

32.  CNAT:  Cjva(M„)  array  of  dimension  NMACH  (1/deg). 

33.  XCPT:  xac(Ma)  array  of  dimension  NMACH  (ft). 

34.  WTAPE:  memg  (0),  empty  fuselage  weight  (lb). 

35.  XCGTP:  x,m  (x-coordinate  of  airframe  (empty)  centre- of- mass 
position,  see  Figure  4,  ft). 

36.  RCGTP:  r«m  (radial  cylindrical  coordinate  of  airframe  (empty) 
centre- of- mass  positions,  see  Figure  4,  ft). 

37.  PCGTP:  <j>„m  (angular  cylindrical  coordinate  of  airframe  (empty) 
centre- of- mass  position,  see  Figure  4,  deg). 

38.  XIXTP:  llem  (slugs. ft2). 

39.  XIYTP:  I",m  (slugs.ft2). 

40.  DIYZ:  Asymmetry  factor  (I?,tm  =  [1  +  DIYZ]  I*,J. 

41.  WPLD:  Weight  payload  (lbs). 

42.  XCGPL:  \PL  (x-coordinate  of  payload  centre- of- mass  position,  see 
Figure  4,  ft). 

43.  RCGPL:  xPL  (radial  cylindrical  coordinate  of  payload  centre- of- mass 
position,  ft). 

44.  PCGPL:  \PL  (angular  cylindrical  coordinate  of  payload  centre-of- 
mass  position,  see  Figure  4,  deg). 

45.  MTYPE:  Number  of  motor  types  (maximum  of  5). 

46.  WTE(I):  Motor  weight  empty  (lb)  of  I- th  motor  type. 

47.  WTPR(I):  Propellant  weight  (lb)  of  I-th  motor  type. 


UNCLASSIFIED 


UNCLASSIFIED 


Al. 


48.  SPIMP(I):  Specific  impulse  of  I-th  motor  type,  IiPi  (sec). 

49.  NTHR(I):  Number  of  thrust  table  points  for  the  I-th  motor  type 
(maximum  of  100). 

50.  CTHR(I):  I-th  motor  type  thrust  correction  factor  (normally  1.0). 

51.  CTHRT(I):  I-th  motor  type  thrust  table  time  correction  factor 
(normally  1.0). 

52.  TTHRT  (J,I):  Time  points  for  I-th  motor  type  (do  not  have  to  be 
equally  spaced,  1  <  J  <  NTHR(I),  secs). 

53.  THRT  (J,I):  Thrust  levels  for  I-th  motor  type  corresponding  to  time 
points  in  TTHRT  (J,I)  (1  <  J  <  NTHR(I),  lbs). 

54.  NKICK:  Number  of  stages  (maximum  of  10). 

55.  NTYPE(I):  Motor  type  used  in  the  I-th  stage. 

56.  NMOT(I):  Number  of  motors  of  type  NTYPE(I)  used  in  the  I-th 
stage,  NMOT(I)  <  3. 

57.  TIG(I):  Time  of  ignition  of  the  I-th  stage  (secs).  The  time  of  ignition 
of  two  successive  stages  may  be  identical,  e.g.  as  would  be  the  case  if 
two  different  motor  types  are  fired  simultaneously. 

58.  XTHR(I):  x-coordinate  (in  F„)  of  the  point  at  which  the  thrust  acts 
for  the  ‘key’  motor  of  the  I-th  stage  (see  Figure  4,  ft).  If  there  is  only 
one  motor  for  the  I-th  stage,  i.e.  NMOT(I)  =  1,  then  the  coordinate 
is  for  that  motor’s  thrust  vector  point  of  action.  If  NMOT(l)  =  2, 
then  the  coordinate  is  for  one  of  the  motor’s  thrust  vectors,  and  the 
program  assumes  that  the  second  motor  is  axisymmetrical  about  the 
x-axis  of  F*  with  respect  to  the  first  motor.  If  NMOT(I)  =  3,  then 
the  coordinate  is  for  one  of  the  motors  off  the  x-axis  of  F„.  The 
program  then  assumes  that  a  second  motor  exists  that  is 
axisymmetrical  about  the  x-axis  of  F„  with  respect  to  the  first  motor, 
and  that  the  third  motor  is  on  the  axis  of  symmetry. 
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59.  RTHR(I):  Radial  cylindrical  coordinate  for  the  I-th  stage  for  the 
same  ‘key’  motor  as  used  for  XTHR(I)  (see  Figure  4,  ft). 

60.  PTHR(I):  Angular  cylindrical  coordinate  for  the  I-th  stage  for  the 
same  ‘key’  motor  as  used  for  XTHR(I)  (see  Figure  4,  deg). 

61.  DTHR(l):  Key  rocket  motor  thrust  vector  reference  frame  Euler 
angle,  see  Appendix  2. 

62.  ATHR(l):  Key  rocket  motor  thrust  vector  reference  frame  Euler 
angle,  see  Appendix  2. 

63.  XEMPT(I):  x-coordinate  (in  F*)  of  the  centre- of- mass  of  the  empty 
motor  case  of  the  ‘key’  motor  of  the  I-th  stage  (see  Figure  4,  ft). 

64.  REMPT(I):  Radial  cylindrical  coordinate  (in  F„)  of  the  centre-of- 
mass  of  the  empty  motor  case  of  the  ‘key’  motor  of  the  I-th  stage  (see 
Figure  4,  ft). 

65.  PEMPT(I):  Angular  cylindrical  coordinate  (in  F*)  of  the  centre-of- 
mass  of  the  empty  motor  case  of  the  ‘key’  motor  of  the  I-th  stage  (see 
Figure  4,  ft). 

66.  XPROP(I),  RPROP(I),  PPROP(I):  Coordinates  analogous  to 
XEMP(I),  REMPT(I),  PEMPT(I)  for  the  cylindrical  coordinates  (in 
Fr)  for  the  centre- of- mass  of  the  propellant  of  the  ‘key’  motor  of  the 
I-th  stage. 

67.  TORQL(I):  Torque  nondimensional  coefficient  about  the  thrust  axis 
of  ‘key’  motor  of  the  I-th  stage  due  to  spiral  grooves  or  vanes  in  the 
nozzei  of  the  motor. 

68.  TCCD:  Coefficient  describing  the  rocket  plume  exhaust  effect  on  the 
drag  of  the  vehicle. 

69.  TCCL:  Coefficient  describing  the  rocket  plume  exhaust  effect  on  the 
aerodynamic  normal  forces  acting  on  the  vehicle. 

70.  TCCP:  Coefficient  describing  the  rocket  plume  exhaust  effect  on  the 
aerodynamic  centre  of  the  vehicle. 
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71.  NALTW:  Number  of  altitudes  at  which  the  wind  velocity  is  given 
(NALTW  <  100).  If  NALTW  =  0,  no  wind  data  is  expected. 

72.  ALTW(I):  Altitudes  above  sea  level  at  which  wind  velocity  vector 
components  are  given  (ft).  These  altitudes  do  not  have  to  be  evenly 
spaced. 

73.  WIXA(I):  W,  corresponding  to  ALTW (I)  (fps). 

74.  WIYA(I):  W2  corresponding  to  ALTW (I)  (fps). 

75.  WIZA(I):  W3  corresponding  to  ALTW(I)  (fps). 

76.  NALTA:  Number  of  altitudes  at  which  atmospheric  temperature  and 
pressure  data  is  given  (NALTA  <  100).  If  NALTA  =  0,  no 
temperature  and  pressure  is  expected. 

77.  ALTA(I):  Altitudes  above  sea  level  at  which  atmospheric 
temperature  and  pressure  data  is  given  (ft).  These  altitudes  do  not 
have  to  be  evenly  spaced. 

78.  TEMPT(I):  The  °7o  deviation  from  standard  temperature 
corresponding  to  ALTA (I). 

79.  PRES (I)  The  °7o  deviation  from  standard  pressure  corresponding  to 
ALTA(I). 

80.  NFIN:  Number  of  cases  (NFIN  <  100). 

A  1.3  SOFTWARE  ALGORITHM 

The  algorithm  used  is  as  follows: 

1)  Input  date  through  unit  NIN  using  subroutine  INPUT. 

2)  Set  NSHOT  =  1.  Total  number  of  cases  for  a  given  run  is  NFIN. 

3)  Using  subroutine  INPUS,  redefine  input  dataset  as  appropriate  for 
case  NSHOT. 
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4)  Output  headers  and  results  as  appropriate  for  value  of  parameter 
NCOUT. 

5)  Compute  starting  parameters,  including  initial  mass  and  inertia 
characteristics  of  the  simulated  vehicle. 

6)  Define  values  of  g,  p,»,  T*  and  g  appropriate  to  lw  with  subroutine 
ATMOS. 

7)  Define  aerodynamics. 

8)  Perform  Runge-Kutta  integration  using  subroutine  RUNK1.  Use 
step  size  DEL  1  if  rocket  motors  are  off,  and  DEL2  if  rocket  motors 
are  on. 

9)  Output  Table  1  data  if  NCOUT  >  0  and  if  an  output  time  point  is 
reached  (at  DELPT  time  intervals). 

10)  Store  results  at  output  time  points  for  all  relevant  parameters  in  the 
array  0(1,  J),  where  I  is  the  parameter  and  J  is  the  output  time  step. 

11)  Compare  outputted  results  with  previous  outputted  results  to  get 
current  maximum  Mach  number,  altitude  and  dynamic  pressure 
data. 

12)  Go  to  Step  6  and  repeat  until  either  —  zf  —  ALTST  -  0.5z,DELT  <  0 
or  t  >  TSTOP.  Here  DELT  =  DELI  or  DEL2  (see  Step  8). 

13)  Output  tabular  data  based  on  data  stored  in  array  0  as  required  by 
parameter  NCOUT.  Convert  to  and  output  in  metric  units  if 
NCOUT  >11,  otherwise  use  Imperial  units. 

14)  If  case  is  finished,  as  determined  by  Step  12,  store  case  summary 
data  (range,  apogee,  flight-time,  apogee-time,  maximum  Mach 
number,  and  maximum  dynamic  pressure).  Output  summary  data 
for  case  just  completed. 

15)  NSHOT  =  NSHOT  +1.  If  NSHOT <  NFIN,  go  to  Step  3. 
Otherwise,  output  case  summary  data  (see  Step  14)  and  STOP. 
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A1.4  GENERAL  NOTES 

A  1.4.1  Integration  Algorithm  Numerical  Instability 

The  Runge-Kutta  algorithm  will  be  unstable  if  DELI  and  DEL2  are  chosen  too 
large,  a  situation  that  will  usually  result  in  an  overflow  condition.  This  problem  may  be 
resolved  by  making  the  integration  step  sizes  DELI  and  DEL2  smaller. 

Al.4.2  Simultaneous  Firing  of  Two  Different  Stages 

The  time  of  ignition  of  two  successive  stages  may  be  identical.  This  would  be 
required  if  two  different  motor  types  are  fired  simultaneously,  or  if  more  than  three 
motors  of  one  type  are  used  in  a  given  stage  (see  the  restrictions  placed  on  NMOT(I)  in 
Section  A  1.2).  Also,  thrust  asymmetries  must  be  handled  by  specifying  the 
characteristics  of  each  motor  in  a  given  stage  individually,  i.e.,  by  treating  them  as 
separate  stage  firings  with  the  same  ignition  time. 

Al.4.3  TCCD,  TCCL,  TCCP 

These  coefficients  are  given  as  inputs  in  order  to  take  into  account  rocket  plume 
exhaust  effects  on  the  aerodynamics  of  the  vehicle  (see  Section  A  1.2).  For 
CRV7/BATS  and  ROBOT-9  type  vehicles,  these  exhaust  effects  on  lift,  side  force  and 
aerodynamic  centre  characteristics  are  negligible,  and  thus  TCCL  =  TCCP  =  1 .0. 

Rocket  motors  fired  in  the  base  of  the  vehicle  will  significantly  reduce  base  drag. 
For  CRV7/BATS  and  ROBOT-9  a  representative  value  for  TCCD  for  such  a  stage  firing 
is  TCCD  =  0.8. 

A  1.4.4  Vehicle  Symmetry  Assumptions 

Whenever  it  is  the  intent  to  use  the  F„'  equations  (NOROL  =1,  see 
Section  A  1.2  and  also  Section  2.2),  all  vehicle  mass,  inertia  and  thrust  characteristics 
must  be  specified  with  axisymmetry  about  the  x-axis  of  F*.  Otherwise,  NOROL  =  0  or 
-1  must  be  used  in  order  to  invoke  the  use  of  the  F*  equations  of  motion. 

For  NOROL  =  0,  axisymmetric  mass,  inertia  and  aerodynamic  characteristics 
must  be  inputted.  General  thrust  characteristics  may  be  used. 
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For  NOROL  =  -1,  no  symmetry  constraints  are  placed  on  inputted  mass, 
inertia,  aerodynamic  and  thrust  characteristics. 

A1.4.S  Output  Tables 

The  output  tables  are  generated  as  determined  by  the  value  of  the  parameter 
NCOUT  (see  Section  A  1.2).  The  contents  of  each  table  may  be  determined  from  the 
output  headers  in  the  subroutine  OUTPU  (see  Section  A  1.5  to  follow). 

Table  1  contains  the  aspect  angle  presentations  of  the  target  relative  to  an 
observer  located  at  a  prespecified  point  relative  to  the  launch  site  reference  frame  F,.  The 
program  assumes  that  the  observer  is  located  at  (xT/,  0,  0)  (see  Section  2.10).  xr,  may 
be  varied  by  changing  the  appropriate  statements  in  the  MAIN  program,  i.e.,  the 
statements  defining  array  elements  0(59,  IOUT)  and  0(60,  IOUT). 
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A1.5  BALSIM  SOURCE  LISTING 

The  following  listing  is  the  IBM  FORTRAN  IV  (H-Extended  Compiler)  version 
of  the  BALSIM  package.  It  is  recommended  that  the  AUTODBL  feature  is  used  so  that 
all  real  arithmetic  is  carried  out  in  double  precision. 
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PROGRAM  B  A  L  S  I  Mi  VERSION  2  (IBM  SYSTEM  370)  OCTOBER  19S2 

COMMON  /INDAT/  NOROL, DELI . TSTOP . ALTST, DELPTi NCOUT, NIOUT. DELS, 
*SLAUN,  SGUID,  SFRLNi  XIO,  YIO,  ALTO.  UO,  VO.  WO.  THETO,  PSIO.  PHIO.  PO  GO  RQ 
*CLP,  CMQ. CLDA. DELTF , B. NMACH, AMACH(20) , CD0T<20), CNAT<20) .  XCPT<  2- 
*WTAPE,  XCGTP.  RCGTP <  PCGTP,  XIXTP,  XIYTP,  DIY2,  WPLD.  XCGPL.  RCGPL.  PCSPu. 
*CLFIN,  CPF  I N.  DDFIN,  PKFIN,  FIFAC,  DXCP.  FXCN.  FXCD, 

*SREF, 

♦MTYPE,  WTE (  5 ) .  WTPR(5),  SPIMPC5),  CTHR<5>.  CTHRT (  5 ) .  NTTHR  (  5  ) . 

*TTHRT( 100,  5). TTHR ( 100,  5),  COR t  <  5 ; - C0R2( 5) , 

»NKICK,  NTYPE  (  10 ) ,  NMGTi  10).  XTHR  (  10) .  ATHR  (  10) ,  DTHR  <  10),  PTHR  (10; 
*XEMPT<  10),  REMPT  (  iO  1  •  PEMPK  10),  XPROP  (  10),  RPROP  (  10),  PPROP  <  10)  ■ 

#T I G ( 10; • RTHR ( 10'  fORQL ( 10 ) • TCCD ( 10 ) , TCCL ( 10) . TCCP< 10) , 

4NALTW,  AlTW(  100),  WIXA(  100),  WIYA1  100) ,  WIZA(  100)  ■ 

■»NALTA,  ALTAI  100),  TMPR  (  100),  PRES<  100), 

»  NFIN 

DIMENSION  INIG( 10 i , CFL ( 10) , CFM< 10  ) , CFN( 10),  CFX  < 10 ) ,  CF Y ( 10 ) ,  CFZ < 10 
*,  WTS  (  10  > .  UITSD<  10) ,  0(60,  200),  A  (3,  3>  -  YPROP  (  10) ,  ZPROP(  10) 

DIMENSION  QBMAX ( 100), XMMAX(IOO), AROSE ( 100), TAPOG(IOO), RANGE ( 100) . 
*FLTIM< 100) , 00(7) 

DATA  CNVRG, CNVGR, GRVO,  NIN, NOUT 
*  /57.  29fe,  0.  017453,  32  174,  5,  6/ 

C 

C  INITIAL  DATA  INPUT 

<  ‘ 

NSHOT  =  0 

CALL  INPUKNIN.  NOUT) 

C 

DO  999  NSHOT  —  1 , NF I N 

C 

c  variable  input  for  EACH  SHOT 

c 

CALL  INPUS (NIN, NOUT, NSHOT, 0, IOUT> 

c 

C  OPTIONAL  ON-LINE  OUTPUT  (NCOUT.  GE.  O) 

c 

IF  ( NCOUT  '  510,511,511 
511  WRITE ( NOUT , 100)  NSHOT 

IF (NCOUT— 10)  520, 521, 521 
521  WRITE (NOUT, 101) 

GOTO  510 

520  WRITE (NOUT, 102) 

510  CONTINUE 

r 

C  INITIAL  CONDITIONS 

C 

U  =  UO  +  0  0001 
V  =  VO 

w  =  wo 

P  »  PO  *  CNVGR 
Q  *  QQ  #  CNVGR 
R  =  RO  *  CNVGR 
PHI  =  PHIC  *  CNVGR 
THET  =  THETO  *  CNVGR 
PS I  =  PSIO  *  CNVGR 
XI  *  XIO 
Y I  =•  YIO 
II  ~  0 

WGTO  -  WTAPE  +  WPLD 
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XCGTO  =  WTAPE*XCGTP  +  WPLD*XCGPL 
YCGTP  -  RCGTP*COS l PCGTP*CNVGR ) 

ZCGTP  =  RCGTP*SIN(PCGTP*CNVGR> 

YCGPL  =  RCGPL*CGS(PCGPL*CNVGR> 

ZCGPL  =  RCGPL*SIN(PCGPL*CNVGR ) 

YCGTQ  =  WTAPE*YCGTP  +  WPLD*YCGPL 
ZCGTD  =  WTAPE*ZCGTP  +  WPLD*ZCGPL 

XIXQ  =  X I XTP*GRVO  +  WTAPE*RCGTP**2  +  WPLD*RCGPL**2 
XIYO  =  X I YTP*GRVO  +  WTAPE* ( XCGTP** 2  +  ZCGTP**2> 

*  +  WPLD*  <  XCGPL**2  +  ZCGPL**2> 

XIZO  =  XIYTPMi  +DIYZ)*GRVO  +  WTAPE* ( XCGTP**2  +  YCGTP**2> 

*  >  WRLD*( XCGPL**2  +  YCGPL**2> 

XIXYG  =- WTAPE* XCGTP* YCGTP  -  WPLD*XCGPL*YCGPL 
XIXZO  =  WTAPE*XCGTP*ZCGTP  +  WPLD*XCGPL*ZCGPL 
XIYZO  =— WTAPE* YCGTP* ZCGTP  -  WPLD* YCGPL*ZCGPL 

C 

DC  500  I  =  1,NKICK 
INIG ( I i  -  1 
DTH=DTHR < 1 )*CNVGP 
SDTH=D3IN(  DTH ) 

CDTH=DC0S ( DTH ) 

ATH  =ATHR< I )*CNVGR 
SATH  =  SIN  ( ATH ) 

CATH  =  COS (ATH) 

PTH  =  PTHR( I )*CNVGR 
SPTH  *  SIN(PTH) 

CPTH  =  COS (PTH) 

CFX(I)  «  1  -FLOAT ( NMOT  < I ) /3 ) /FLOAT ( NMOT ( I ) ) 

CFX(I)  =  CDTH*CATH*CFX ( I )+l.  -CFX ( I ) 

CFY  (  I )  ~-CDTH*SATH*SPTH-SDTH*CPTH 
CFZ ( I )  =-CDTH*SATH*CPTH+SDTH*SPTH 
CFL(I)  =  TORQL ( I ) *CFX  <  I  ) 

CFM(I)  =  CFZ ( I ) *XTHR ( I  )  -  RTHR ( I ) *SPTH*CFX ( I )  *  TORQL ( I ) *CFY ( I) 
CFN'I)  =— CFY ( I ) *XTHR ( I )  -  RTHR ( I ) *CPTH*CFX ( I )  +  TORQL ( I ) *CFZ ( I ) 

r* 

U 

K  =•  NTYPE(I) 

FNMOT  =  FLOAT ( NMOT  < I )  ) 

F3M0T  =  FLOAT ( 1 +NMQT  < I > /2 ) /FNMOT 
YPROP(I)  =  RPROP ( I >  *CQS ( PPROP ( I >  *CNVGR ) 

ZPROP ( I )  =  RPROP ( I )*SIN(PPRQP ( I ) *CNVGR ) 

YEMPT  -  REMPT ( I >  *COS ( PEMPT ( I ) *CNVGR ■ 

ZEMPT  =  REMPT'.  I  >*SIN(PEMPT(  I  }*CNVGR  :■ 

WTS  < I )  -  WTPR ( K ) *FNMOT 
WTSD  (  I  >=0.  0 

WGTO  =  WGTO  +  WTE(K)*FNMOT 
XCGTO  »  XCGTO  *  WTE < K > *XEMPT ( I > *FNMOT 
XIXO  *  XIXO  +  WTE • ft ) *REMPT ( I ) **2*FNM0T *FSMOT 
XIYO  =  XIYO  *  WTE ( K )  *  t  XEMPT ( I > **2  *  ZEMPT**2*FSM0T ) *FNMOT 
XIZO  =  XIZO  +  WTE(K)*(  XEMPT'.  I  )**2  +  YEMPT*  *2*FSM0T  >  *FNMOT 
IF(NMOT(  I  ,-i  >  519.519,500 
519  CONTINUE 

YCGTO  =  YCGTO  +  WTE ( K ) *YEMPT*FNMOT 
ZCGTO  =  ZCGTO  +  WTE ( K )* ZEMPT  *FNMOT 
XIXYO  =  XIXYQ  -  WTE ( ft  > * XEMPT ( I ) *YEMPT*FNMOT 
XIXZO  =  XIXZO  +  WTE ( ft  > *XEMPT ( I >  *ZEMPT*FNMOT 
XIYZO  *  XIYZO  -  WTE(K)*YEMPT*ZEMPT*FNMOT 
500  CONTINUE 

INITIAL  AERODYNAMICS 

PHIFR  --  PHFIN*CNVGR 
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SIPHF  -  SINK  PHI FR , 

COP HR  -  COS ( PHIFR ; 

CYFIO  --  -CL.FIN+DDFIN+SIPHF 
CZFIO  =  -CLFIN*DDFIN*CQPHF 
CYFIA  =  -FIFAC*CLFIN*SIPHF*CNVRG 
CZFIA  =  -FIFAC*CLFIN*COPHF*CNVRG 


KUTTA 


RHFCT  =  I 
VSFCT  =  1 

IF (NALTA-1 )  501;  502;  504 
CONTINUE 
GOTO  504 
N1  *  O 

CALL  ATMOSlRHO.  PP,  TK,  VS,  GRAV,  ALTO,  TMPR,  PRES,  ALTA,  Nl,  RHFCT ■  VSFCT.’ 
RHFCT  =  PRESd  >*TK/(PP*TMPR(  l  >  > 

VSFC T  =  SORT  <  TMPR ( 1 ) /TK ) 

CONTINUE 

GBMAX ( NSHOT )  =  Q 

XMMAX( NSHOT)  =  0. 

APOGE ( NSHOT )  =  0 
TIME  =  O. 

PTIME  =  0. 

XI ALP  =  0 
iour  =  O 
KILL  *  0 
DELT  =  OElO 
KUTTA  =  4 
GOTO  4 

KUTTA  =  KUTTA  +  1 
GOTO  (3,  2,  3,  2)  •  KUTTA 
CONTINUE 

TIME  *  TIME  +  DELT  *  0.  5 
CALL  RUNK1  (KUTTA.-  OELT ,  U,  UD,  1  ) 

CALL  RUNK1 (KUTTA, DELT,  V,  VD,  2) 

CALL  RUNKI (KUTTA,  DELT,  W,  WD.  3) 

CALL  RUNK1 (KUTTA,  DELT, P, PD,  4) 

CALL  RUNK1 (KUTTA,  DELT,  G,  GD,  5) 

CALL  RUNKi  (KUTTA,  DELT,  R,  RD,  6) 

CALL  RUNKI (KUTTA.  DELT,  PHI,  PHID-  7) 

CALL  RUNKI (KUTTA, DELT, THET,  THETD,  8) 

CALL  RUNKI  (KUTTA,  DELT,  PSI.  PSID,  9? 

CALL  RUNKI  (KUTTA,  DELT,  XI,  XDI,  10) 

CALL  RUNKI (KUTTA, DELT,  Y I,  YD  I,  11 ) 

CALL  RUNKI (KUTTA, DELT,  Z I,  ZD I.  12) 

DO  13  I  *  l.NKICK 
J  =  1+12 

CALL  RUNKI  (KUTTA, DELT,  WTS< I ) ,  WTSD( I > ,  J) 

CALL  RUNKI  (KUTTA,  DELT,  XIALP,  XIALPD, 24) 

CONTINUE 

PHI  =  AMOD(PHI, 6  28318531 ) 

STH»S I N ( THET ) 

XIALPD  =  GRVQ*STH 
CTH*COS(THET> 

SPS=SIN  <  PSI ) 

CPS«COS(PSI> 

SPH«SIN(PHI ) 

CPH«COS(PHI> 

A(l. 1)«CTH*CPS 
Ad,  2>«CTH*SPS 
Ad,  3>»-STH 

A <2, 1 )*SPH*STH*CPS-SPS*CPH 


CALL 

CALL 

CALL 

CALL 

CALL 
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A<2, 2)=SPS*STH*SPH+CPS*CPH 
A (2, 3>=SPH*CTH 
A(3,  1 >=CPS*CPH*STH+SPS*SPH 
A  ( 3,  2 ) =SPS*CPH*STH-CPS*SPH 
A (3. 3)=CPH*CTH 

XDI  =  A( 1. 1 >*U+A<2, 1 )*V+A(3, 1 >*U 
YDI  =  A< 1, 2»*U+Al2, 2>*V+A(3, 2)*W 
ZDI  =  A<  1, 3>*U+A<2,  3>*V+A(3,  3>*W 

ALT=-ZI  +  ALTO 

CALL  ATMOS (RHO,  POMP, TKELV,  VS, GRAV, ALT. TMPR, PRES, ALTA,  NALTA. 
♦RHFCT. VSFCT) 

URL  =  U 
VRL  *  V 
URL  *  W 

IF<NALTU) 14, 18, 14 

14  CALL  WINDtWIXA,  WIYA,  WIZA,  ALTW,  NALTW,  ALT,  A,  WX,  WY,  WZ  ) 

URL  =  U-WX 
VRL  =  V-WY 
URL  -  U-UZ 

IS  VR  =  SQRT<URL*URL+VRL*VRL+URL*URl^ 

XMACH=VR/VS 
GBAR=0  5»RH0*VR*VR 
CDRAG  =  10 
IF  <  I  PASS  >  913,913,912 

913  ALPHA=  O  0 
BETA  =  0.  0 
ALPHD  =  0.0 
BET AD  =00 
IPASS  =  1 
GO  TO  914 

912  ALPHA  =  ATAN  (URL/URL) 

BETA  =  ATAN  ( VRL /SORT < ORL*URL+WRL*URL > > 

ALPHD= ( UD/VR-UD*UD/ <  VR*VR ) ) *COS ( ALPHA  >  *COS ( ALPHA ) 

TEMP  =  SORT  <URL*URL+URL*WRL> 

BETAD= ( TEMP*VD-VRl*  <  URL*UD+URL*WD ) / TEMP ) / ( TEMP*TEMP*COS ( BETA , i 
UF IN  =  WRL*COPHF  +  VRL*SIPHF 
ALFIN  =  ATAN (UF IN/URL) 

914  CONTINUE 

C 

C  THRUST, MASS, C  G, MOM. IN  ARE  COMPUTED 

C 

GOTO  <  390, 300, 390, 300 ) , KUTTA 
300  FEX  =  O 
FEY  =  O 
FEZ  =  O. 

EML  =  0. 

EMM  =  O. 

EMN  =  0. 

LOMO  =  0 
THRT  =  0. 

XIDX  =  0. 

XIDY  =  O. 

XIDZ  =  0. 

SWD  =  0. 

SWDX  =  0 
SWDY  =  0. 

SWDZ  =  0. 

8UDX2  =  0. 

8UDY2  =  0. 

8WDZ2  =  0. 

LONSY  =  0 
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IF  (NOROL )  310,311,311 

310  CONTINUE 
XIDXY  *  0. 

XIDXZ  «  0. 

XIDYZ  *  0. 

SWDXY  *  0. 

SWDXZ  =  0. 

SWDYZ  =  0. 

3 11  CONTINUE 
C 

WGHT  a  WGTO 
XCGT  *  XCGTO 
XIX  *  XIXO 
XIY  *  XIYO 
XIZ  =  XIZO 

YCGT  =*  YCGTO 
ZCGT  a  ZCGTO 
X IXY  a  XIXYO 
XIXZ  *  XIXZO 
XI YZ  a  XIYZO 

SUM  OVER  ALL  KICKS 

DO  308  Iai,NKICK 
WTSD( I )»0.  0 

IF<TIME-TIG(I> >  301,302,302 
302  K=NTYPE(I) 

LaNTTHR(K) 

IF<TIME~TIG< I l-TTHRT <L,  K  >*CTHRT ( K > >304,  304,  301 
304  T  a(TIME-TIG< I ) ) /CTHRT (K) 

LOMO  *  1 
LONSY  =  2 

FSMOT  *  FLOAT  < 1 +NMOT  < I ) /2 ) /FLOAT  (  NMOT ( I ) ) 

CALL  INTPL  (TTHR ( 1,  K) ,  TTHRT < 1,  K) ,  THR,  T,  INIG ( I > ,  L  > 
THR  =  THR  *  FLOAT ( NMOT ( I ) ) #CTHR ( K ) 

THRT  =  THRT  +THR 

WTSD(I)  *  -THR/ (CTHRT ( K ) #SP IMP ( K )*CTHR (K ) ) 

XMASD  a  -WTSD< I J/GRVO 
SMASD  »  XMASD*FSMOT 
CDRAG  a  TCCD(I) 

FEX  *  FEX  +  THR#CFX  < I ) 

EML  *  EML  +  THR*CFL( I ) 

SWD  «  SWD  WTSD(I) 

DWDX  *  WTSD ( I ) *XPROP  < I > 

DWDY  «  WTSD ( I ) *YPROP ( I > 

DWDZ  «  WTSD  < I >  *ZPROP ( I > 

8WDX  «  SWDX  +  DWDX 
SWDY  •  SWDY  ♦  DWDY*FSMOT 
SWDZ  a  SWDZ  ♦  DWDZ*FSMOT 
SWDX2  >  SWDX2  DWDX*XPROP ( I ) 

SWDY2  »  SWDY2  +  DWDY*YPROP (I > *FSMOT 
SWDZ2  =  SWDZ2  +  DWDZ#ZPROP < I >*FSMOT 
IF(NMOT ( 1 )  —  1  )  301,307.301 
307  FEY  a  FEY  +THR*CFY ( I > 

FEZ  ■  FEZ  +THR*CFZ  < I ) 

EMM  =  EMM  +THR*CFM< I > 

EMN  =  EMN  +THR*CFN(I) 
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SWDXY  =  SWDXY  +  DWDX*YPROP < I ) 

SWDXZ  =  SWDXZ  +  DWDX*ZPRQP(I) 

SWDYZ  =  SWDYZ  +  DWDY*ZPROP < I > 

LONSY  *  1 
301  CONTINUE 

WGHT  =  WGHT  +  WTS ( I ) 

XCGT  =  XCGT  +  WTS ( I ) *XPROP ( I ) 

XIX  =  XIX  +  WTS ( I ) *RPRQP  < I > **2*FSM0T 

X I Y  =  X I Y  +  WTS< I )*( XPROP ( I )**2  +  ZPROP ( I > **2*FSM0T > 

XIZ  =  XIZ  +  WTS  (  I  )  *  <  XPROP  (  I )  **2  +  YPROP < I >**2*FSM0T> 

IF  < NMOT  V  I  >  —  1  >  320,320,308 
320  X I X Y  a  XIXY  -  WTS (I)*  XPROP ( I ) *YPROP  < I ) 

XIXZ  =  XIXZ  +■  WTS  C  I )  *  XPROP  (  I  )*  ZPROP  ( I  ) 

XIYZ  =  XIYZ  -  WTS <  I  )*  YPROP  ( I ) «•  ZPROP (  I  ) 

YCG'T  *  YCGT  +  WTS <1 > *YPROP <1 > 

ZCGT  =  ZCGT  +  WTS< I >*ZPRGP< I > 

30S  CONTINUE 

IF  ( NORQL  )  331,330,330 

SYMMETRICAL  BODY 

330  XMASS  =  WGHT/GRVG 
YCGT  =  O. 

ZCGT  *  O 

XCGT  =  XCGT /WGHT 

XIX  =  XIX/GRVO 

XI Y  »  XIY/GRVQ  -  XMASS* XCGT**2 
XIZ  =  XIZ/GRVG  -  XMASS*XCGT**2 
X IDX  a  SW0Y2  +  SWDZ2 

X IDY  =  SWDX2  +  SWDZ2  +  SWD*XCGT**2  -  2.  *SWDX*XCGT 

X IDZ  =*  SWDX2  +  SWDY2  +  SWD*XCGT**2  -  2.  *SWDX*XCGT 

X IDX  =  X IDX/GRVQ 

X IDY  =  XIDY/GRVO 

XIDZ  =  XIDZ/GRVO 

•  EMM  =  EMM  -  FEZ*XCGT 
EMN  =  EMN  +  FEY*XCGT 
GOTO  390 


UNSYMMETRICAL  BODY 

331  XMASS  =  WGHT/GRVO 
XCGT  *  XCGT /WGHT 
YCGT  =  YCGT /WGHT 
ZCGT  =  ZCGT /WGHT 
CGZ  =  XCGT**2  *  YC,GT**2 
CGY  =  XCGT  **2  +  ZCGT  **2 
CGX  *  YCGT#*2  +  ZCGT**2 
CGXY  =-  XCGT*YCGT 
CGXZ  a  XCGT*ZCGT 
CGY7  =  -YCGT*ZCGT 

XIX  =  XIX/GRVO  -XMASS*CGX 
XI Y  =  X  I Y/GRVO  -XMASS*CGY 
XIZ  =  X I Z/GRVO  -XMASS*CGZ 
XIXY  =  +  X I X Y/GRVO  -  XMASS*CGXY 
XIXZ  =  +  X I XZ/GRVO  -  XMASS*CGXZ 
XIYZ  =  +  X IYZ/GRVO  -  XMASS*CGYZ 
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SWCX  =  SWDX*XCGT 
SWCY  =  SWDY*YCGT 
SWCZ  =  SWDZ*ZCGT 

XIDX  *  SWDY2  +  SWDZ2  ♦  SWD*CGX  -  2.  *  <  SWCY+SWCZ  > 

XIDY  =  SWDZ2  +  SWDY2  +  SWD*CGY  -  2.  *  <  SWCX+SWCZ  > 

X IDZ  =  SWDY2  +  SWDZ2  +  SWD*CGZ  -  2.  *<SWCX+SWCY> 

XIDXY  —SWDXY  +  SWD*CGXY  +  XCGT*SWDY  +  YCGT*SWDX 

XIDXZ  »  SWDXZ  +  SWD*CGXZ  -  XCGT*SWDZ  -  ZCGT*SWDX 

XIDYZ  —SWDYZ  +  SWD*CGYZ  +  YCGT*SWDZ  +  ZCGT*SWDY 

XIDX  =  XIDX/GRVO 
XIDY  «  XIDY/GRVQ 
XIDZ  a  XIDZ/GRVO 
XIDXY  »  X IDXY/GRVO 
XIDXZ  *  XIDXZ/GRVC/ 

XIDYZ  =  X IDYZ/GRVO 

CALCULATION  OF  INVERSE  TENSOR  OF  INERTIA 

YIX  *  XIY*XIZ  -  XIYZ**2 
YIXY  »  XIXZ*XIYZ  +  X I XY*X  I Z 
YIXZ  *  X I XY*X I YZ  +  X I  XZ*X I Y 
YNEN  =  X I X*YI X  +  XIXY*YIXY  -  XIXZ*YIXZ 
YIX  a  YI X/YNEN 
YIXY  =  YIXY/YNEN 
YIXZ  =  YIXZ /YNEN 
YIY  a  <XIX*XIZ  -  XIXZ**2)/YNEN 
YI YZ  =  <XIXZ*XIXY  +  X I YZ*XIX ) /YNEN 
YIZ  a  (XI X*X I Y  -  X I X  Y**2 ) /YNEN 

THRUST  MOMENTS 

EML  *  EML  -  FEZ*YCGT  -  FEY#ZCGT 

EMM  *  EMM  -  FEZ*XCGT  +  FEX*ZCGT 

EMN  a  £MN  +  FEY*XCGT  +  FEX*YCGT 

390  CONTINUE 

AERODYNAMICS 

GSREF  a  QBAR*SREF 
GSVRF  *  QSREF*B*B*CNVRG/(2.  *VR  > 

CALL  INTEQ  (CDOT,  CNAT,  XCPT,  AMACH,  CDO,  CNA,  XCP,  XMACH,  NMACH > 

C 

XCP  -XCP  +  DXCP 

FAX  =  -CDO  *QSREF*CDRAG*FXCD 

FAY  a  -CNA  *BETA*QSREF*CNVRG*FXCN 

FAZ  *  -CNA*ALPHA*QSREF*CNVRG*FXCN 

DDF AY  »  QSREF* ( C YF I A*ALF I N+C YF I 0 ) 

DDFAZ  =  QSREF  * ( C  ZF I A*ALF I N+C  ZF 1 0  > 

C 

AL  a  CLP*P*GSVRF  ♦  CLDA*DELTF*QSREF*B  -  FAY*ZCGT  -  FAZ*YCGT 
AM  a  FAZ* ( XCP-XCGT )  +  CMQ*Q*QSVRF  +FAX*ZCGT  +  DDFAZ*<CPFIN-XCGT> 
AN  «»FAY*< XCP-XCGT)  +  CMQ*R*QSVRF  +FAX*YCGT  -  DDFAY* (CPF IN-XCGT ) 
C 

FAY  -  FAY  +  DDFAY 
FAZ  -  FAZ  *  DDFAZ 
C 

FXT»FAX+FEX 

FYT»FAY*F£Y 
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FZT*FAZ+FEZ 
TL  =  AL  +  EML 
TM  *  AM  +  EMM 
TN  =  AN  +  EMN 


DYNAMICS 

THETD  =  G*CPH-R*SPH 
PSIDa ( Q*SPH+R*CPH ) /CTH 
I F  <  NOROL )  403. 40 1 .  402 

401  PHID  *  P  *  PSID*STH 
UD=FXT/XMASS-W*Q+V*R+GRAV*A( 1,  3) 
VD»FYT/XMASS+W#P-U*R+GRAV*A(2,  3) 
UD»FZT/XMASS+U*Q-V*P+GRAV*A (3.3) 

HCA  *  <XIZ-XIY)*Q*R 
HCB  =  (XI X-X I Z )*P*R 
HCC  *  <  XIY-XIX  )*P*G 

HDA  »  X IDX*P 
HDB  »  X IDY*Q 
HDC  a  XIDZ*R 

PD  =  (TL-  HDA  -  HCA) /XIX 
GD  »  (TM-  HDB  -  HCB  > /X I Y 
RD  a  (TN-  HDC  -  HCO/XIZ 
GOTO  410 

402  PHID  »  PSID*STH 

UD  a  FXT/XMASS-W*Q+V*R+GRAV*A(1, 3) 


VD  *  FYT/XMASS  - 
UD  a  FZT/XMASS  ♦ 

HCA  =  O. 

HCB  =  0 
HCC  »  O. 

HDA  a  P#X IDX 
HDB  »  Q*X IDY 
HDC  »  R*XIDY 

PD  *  <Tl.-  HCA  -  ) 
GD  a  <TM-  HCB  -  » 
RD  ■  (TN-  HCC  -  » 
GOTO  410 


0RAV*A(2. 3) 
GRAV*A ( 3.  3) 


HDA ) / X I  X 
HDB  > / X I Y 
HDC ) /XI Y 


403  PHID  ■  P  +  PSID*STH 

UD  «  FXT/XMASS  -W*Q  ♦  V*R  ♦  GRAV*A(1,3> 
VD  »  FYT/XMASS  +W*P  -  U*R  +  GRAV*A(2.3) 
UD  «  FZT/XMASS  U*G  -  V*P  ♦  0RAV*A(3.3> 

HA  »  P*XIX  -  G*XIXY  -  R*XIXZ 
HB  *-P#XIXY  ♦  Q#XIY  -  R#XI YZ 
HC  «-P*XIXZ  -  G*X I YZ  ♦  R*XIZ 

HCA  >  G*HC  -  R*HB 
HCB  -  R*HA  -  P*HC 
HCC  ■  P*HB  -  Q*HA 

HDA  **  P*XIDX  -  Q*XIDXY  -  R*XIDXZ 
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HDD  =-P*XIDXY  +  G*X IDY  -  R*XIDYZ 
HDC  =-P*XIDXZ  -  Q*X IDYZ  +  R*XIDZ 
C 

QMA  =  TL  -  HDA  -  HCA 

OMB  =  TM  -  HDB  -  HCB 

OMC  =  TN  -  HDC  -  HCC 

C 

PD  =  OMA*YIX  +  OMB*YIXY  +  OMC*YIXZ 

GD  =  OMA*YIXY  +  OMB*YIY  +  OMC*YIYZ 

RD  =  OMA*YIXZ  +  OMC*YIYZ  +  OMC*YIZ 

GOTO  410 
410  CONTINUE 
C 

C  KINEMATIC  RESTRICTIONS  WHILE  BIRD  IS  ON  LAUNCHER 

S  =  SORT ( ( X I-XI0>**2+<  YI-YI0)*#2+ ( ZI )**2> 

IF (S-SGUID )  202.202,230 
202  GD  =  00 
RD  *  0.  0 
PD  -  0.  0 
V  «  VO 
W  =  WO 
VD  —  0.  O 

UD  =  AMAX1  (UD,  0.  ) 

WD  =  O.  O 
230  CONTINUE 

IF  (KUTTA— 4)  1,5,5 
5  KUTTA  =  O 
C 

C  CALCULATION  OF  TRAJECTORY-PARAMETERS 

C 

IF (GBAR-QBMAX (NSHOT ) )  255, 255, 254 
254  QBMAXt NSHOT)  =  GBAR 

255  IF (XMACH-XMMAX< NSHOT ) )  257,257,256 

256  XMMAXt NSHOT)  »  XMACH 

257  IF ( ALT- APOGE< NSHOT )  )  259.259,258 

258  APOGE( NSHOT)  ■  ALT 
TAPOG( NSHOT)  ■  TIME 

259  CONTINUE 
C 

C  CALCULATION  OF  DELT 

C 

IF(LOMO)  231.231,232 

231  DELT  *  DELI 
GOTO  233 

232  DELT  *  DEL2 

233  CONTINUE 
C 

C  END  OF  RUNGE-KUTTA  LOOP 

C  END  OF  SHOT? 

C 

IF(ZDI)  240,240,241 

241  IF<-ZI-ALTST-ZDI*DELT*0.  5)  290,290,240 
240  IF(TIME-TST0P+DELT*0.  5)  250,290,290 
C 

C  OUTPUT  WANTED'* 

C 

250  IF< TIME— PTIME+.  001*DELT)  1,251,251 

251  I OUT  »  IOUT+1 
PTIME  -  PTIME+DELPT 

IF(  I  OUT- 199)  252,290,290 
C 
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C  OUTPUT  ARRAY  IS  FILLED 

C 

252  0(1, I OUT )  =  TIME 
0(2, I OUT )  =  XI 
0(3, I OUT )  *  YI 
0(4.  I OUT )  =  -  ZI 

C 

C  -ZI  =  ALTITUDE  CENTRE-OF-MASS  ABOVE  GROUND  LEVEL. 

C 

0(5,  I OUT )  =  SORT  <U**2+V**2+W**2) 

TEMP  =  SORT (XD 1**2+ YD 1**2) 

0(6, IOUT)  =  CNVRG  *  ATAN2 < -ZDI . TEMP > 

0(7, IOUT)  =  CNVRG  *  ATAN2 ( YDI , XDI ) 

0(8, IOUT)  =  UD 
0(9, IOUT)  *  VD 
0( 10, IOUT)  =  WD 
0(11,  IOUT)  =  THET*CNVRG 
0(12. IOUT)  =  PSI  *CNVRG 
0(13, IOUT)  «  PHI  *CNVRG 
0(14. IOUT)  ■  P  *CNVRG 
0(15.  IOUT)  =  Q  +CNVRG 
0( 16,  IOUT)  =  R  *CNVRG 
0(17, IOUT)  =  PD  *  CNVRG 
0(18, IOUT)  *  QD  *CNVRG 
0(19, IOUT)  »  RD  *CNVRG 
0(20,  IOUT)  = ALP HA*C NVR G 
0(21, IOUT)  =  BETA*CNVRG 
0(22, IOUT)  «  QBAR 
0(23, IOUT)  =  XMACH 
0(24. IOUT)  =  FAX 
0(25, IOUT)  =  FAY 
0(26,  IOUT)  «  FAZ 
0(27. IOUT)  *  FEX 
0(28, IOUT)  -  FEY 
0(29, IOUT)  *  FEZ 
0(30,  IOUT)  =  AL 
0(31, IOUT)  =  AM 
0(32, IOUT)  =  AN 
0(33,  IOUT)  *  EML 
0(34, IOUT)  =  EMM 
0(35, IOUT)  *  EMN 
0(36, IOUT)  =  WGHT 
0(37, IOUT)  =  XCGT 
0(38, IOUT)  *  XCP 
0(39. IOUT)  =  XIY 
0(40,  IOUT)  =  XIX 
0(41, IOUT)  =  YCGT 
0(42, IOUT)  =  ZCGT 
0(43, IOUT)  =  XIZ 
0(44, IOUT)  =  XIXY 
0(45, IOUT)  =  XIXZ 
0(46, IOUT)  =  XIYZ 
0(47,  IOUT)  =»-HCA 
0(48, IOUT)  =-HCB 
0(49, IOUT)  =-HCC 
0(50, IOUT)  =-HDA 
0(51, IOUT)  «-HDB 
0(52, IOUT)  *-HDC 
0(53, IOUT)  *  WX 
0(54, IOUT)  »  WY 
0(55, IOUT)  *  WZ 
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WW1=WX*A( 1, 1 )+WY*A<2, 1 )+WZ*A<3, 1 ) 

WW2=WX*A ( 1 ,  2)+WY*A(2, 2)+WZ*A<3,  2) 

WW3=WX*A ( 1 , 3>+WY*A(2,  3>+WZ*A(3,  3> 

0(56. I OUT )  =  WW1 
0(57,  I OUT  >  *  WW2 
0(58. I OUT)  -  WW3 

0(59,  I0UT)  =  CNVRG*ATAN(-ZI/ ( 18500.  /.  3048-XI ) ) 

0(60,  I  OUT )  -  CNVRG*ATAN(-YI  /  ( 18500.  /  3048-XI )) 

OPTIONAL  OUTPUT  ON-LINE  FOR  NCOUT.  GE.  0 

IF  ( NCOUT  >  269,261,261 

261  IF  (NCOUT— 10)  263,262,262 

262  00(1 )  -  XI*. 3048 

00(2)  =  -ZI*.  3048 
00(3)  =  0(5, IOUT)*. 3048 
00(4)  *  QBAR*47.  9 
00(5)  *  FEX*4.  45 
00(6)  -  FAX*4.  45 
00(7)  =  WGHT*.  453 

WR ITE( NOUT ,  110)  TIME,  (00(1),  1  =  1, 3), 0(6,  IOUT) , 0(20-  IOUT) , 0(14,  IOUT) 
*,  XMACH,  (00(1),  1=4,  7) 

GOTO  269 

263  WRITE ( NOUT ,  1 10)  TIME,  XI-  ALT,  0(5,  IOUT),  0(6-  IOUT),  0(20,  IOUT), 

*0(14,  IOUT),  XMACH,  Q8AR,  FEX,  FAX,  WGHT 

269  CONTINUE 

IF  (KILL)  1,1,299 
290  KILL  =  1 

IOUT  =  IOUT  *  1 
GOTO  252 

END  OF  SHOT 

299  ICOUT  *  NCOUT 

CALL  OUTPU ( 0, IOUT. NOUT, ICOUT, NSHOT) 

RANGE (NSHOT)  -  XI 
FLTIM( NSHOT)  ■  TIME 

CALL  OUT AB ( NOUT , NCOUT. NSHOT, RANGE, FLTIM, APOGE, TAPOG, QBMAX, XMMAX,  1 ) 
999  CONTINUE 

END  OF  RUN 

CALL  OUT  AB ( NOUT ,  NCOUT, NFIN  , RANGE, FLTIM, APOGE, TAPOG.  QBMAX,  XMMAX,  2) 
STOP 

100  FORMAT ( 1 H 1 ,  X ,  40 (  1 H* ) ,  '  RESULTS  SHOT  NO  ',  13.  X.  40<  1H* ) ,  //, 

*3X.  'TIME  ' ,  6X,  'RANGE',  X,  '  ALTITUDE X,  VELOCITY ',  4X,  'ELEV',4X. 

*  'ALPHA ' ,  X,  'ROLLRATE',  5X,  'MACH'.X,  'DYN.  PRES',3X,  'THRUST',  5X.  'DRAG', 
*4X»  'WEIGHT') 

101  FORMATOX,  '  (SEC)  ',  4X.  '  (M)  ',  5X,  '  (M)  ',  5X,  '(M/S)  '.  5X,  '(DEG)  ',  4X. 

* '  ( DEG )  ',  2X,  '(DEO/S)  ',  5X,  '(  )  ',  3X,  '  (N/M2)  ',  5X.  '  (N)  '.  7X,  '  (N)  ',  5X, 

*  '  (  KO )  ',  /) 

102  FORMATOX,  '  (SEC)  ',  3X,  '  (FT)  ',  4X,  '  (FT)  ',  4X,  '  (FT/S)  ',  5X,  '  (DEG)  ',  4X, 

*  '  ( DEG )  ',  2X,  '(DEG/S)  ',  5X,  '(  )  ',  2X,  '  (LB/FT2)  ',  4X,  '  (LB)  ',  6X,  '  (LB)  ', 
*4X,  '(LB)  ',  /) 

110  FORMAT (F8.  2,  F9.  0.  F9.  0,  F9.  1,  F9.  2,  F9  2,  F9.  0,  F9.  2,  F9.  0,  F9  1.F9  1,  F9  2 
*) 

END 
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SUBROUTINE  INPUS(NIN, NOUT, NSHOT ,  O, IOUT) 

COMMON  /INDAT/  NOROL,  DEL  1 ,  TSTOP*  ALTST , DELPT.  NCOUT <  NIOUT , DEL2 
♦SLAUN,  SGUID,  SFRLN.  XlO.  YIO.  ALTO,  UO,  VO,  WO,  THETO,  PSIO,  PHIO,  PO,  QO,  RO, 
♦CLP,  CMG> CLDA,  DELTF,  B,  NMACH, AMACH < 20 ) , CDOT ( 20 ) , CNAT < 20  > , XCPT<20), 
♦WTAPE. XCGTP, RCGTP, PCOTP, XIXTP, XIYTP, DIYZ. WPLD, XCGPL,  RCGPL, PCGPL. 
♦CLFIN,  CPFIN,  DDF  IN,  PHFIN,  FIFAC.  DXCP,  FXCN,  FXCD. 

♦SREF, 

♦MTYPE, WTE (5 ) , WTPR  <  5 ) , SP IMP ( 5 ) , CTHR ( 5) , CTHRT  <  5  >  >  NTTHR ( 5  > , 

♦TTHRT <  100,  5),  TTHR  ( 100,  3) ,  C0R1  (  5) ,  C0R2(3> , 

♦NKICK, NTYPE( 10), NMOTI 10), XTHR ( 10 ) , ATHR ( 10) , DTHR <10 ) ,  PTHR ( 10 ) . 
♦XEMPT< 10), REMPT(IO), PEMPT(IO), XPROP < 10 > ,  RPROP < 10 ) ,  PPROP ( 10) , 

♦TIG  < 10), RTHR(IO), TORGL< 10), TCCD( 10 ) ,  TCCL < 10) ,  TCCP < 10) , 

♦NALTW, ALTW(IOO),  WIXA(IOO),  WIYA(IOO).  WIZA(IOO), 

♦NALTA,  ALTAI  100) ,  TMPR< 100),  PRES < 100), 

♦  NFIN 

DIMENSION  0(60,200) 

IF  (NSHOT.  EG  1)  THET0*45.  0 
IF  (NSHOT.  EG.  2)  THET0=3C  0 
IF (NSHOT  EG  3)  THET0=55.  0 
IF  (NSHOT.  EG.  4)  TH£T0*60.  0 
IF  (NSHOT.  EG.  3)  THET0*63.  0 
IF  (NSHOT.  EG.  6)  THET0=70. 0 
WRITE (6, 1 )  NSHOT, THETO 

F0RMAT(1H1,  'CASE*  ',  13,  'THETQ(DEG)  =  '.1PG13.  3) 

RETURN 

END 
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SUBROUTINE  INPUT  (NIN. NOUT) 

COMMON  /INDAT/  NORQL. DELI. TSTOP, ALTST, DELPT, NCQUT, NIOUT, DEL2. 
#SLAUN.  SGUID,  SFRLN.  XIO,  YIO,  ALTO-  UO,  VO,  WO,  THETO,  PSIO.  PHIO,  PO.  QO,  RO. 
♦CLP, CMQ, CLDA, DELTF, B, NMACH, AMACH ( SO ) ,  COOT ( 20  > ,  CNAT<20>  .  XCPT(2G), 
♦WTAPE, XCGTP, RCGTP, PCGTP, XIXTP, XIYTP. 01 YZ, WPLD, XCGPL, RCGPL- PCGPL, 
♦CLFIN, CPFIN, DDF IN, PHFIN, FIFAC, DXCP,  FXCN, FXCD. 

♦SREF, 

♦MTYPE,  WTE  ( 3  > ,  WTPR  (  5 ) ,  SPIMP  (3),  CTHRO),  CTHRT  (  3 ) ,  NTTHR  <  5  ) , 

♦TTHRTt  100,  3),  TTHR  < 100,  3),  CORKS) .  C0R2<  3), 

♦NKICK,  NTYPE ( 10 ) ,  NMOT  ( 1 0 ) .  X THR  (10),  ATHR  <  1 0 ) ,  DTHR  <  1 0  > ,  PTHR  (  10  > . 
EXEMPT < 10) , REMPT( 10) , PEMPT< 10),  XPRQP ( iO > . RPROP ( 10 ) ,  PPROP < 10 ) , 

♦TIG (  10),  RTHR<10),  TORQL(IO),  TCCD(  iO),  TCCL<  10),  TCCPUO), 

♦NALTW,  AlTW<  100),  WIXA(  100),  WIYA(  100) ,  WIZAUOO). 

♦NALTA, ALTA< 100),  TMPR ( 100),  PRES< 100) , 

♦  NFIN 


READ(NIN, 100) 
READ (NIN, 101 ) 
READ (NIN,  101 ) 
READ (NIN,  101 ) 
READ (NIN,  101 ) 
READ (NIN.  101 ) 
READCNIN,  101) 
READ (NIN,  101 ) 
READ (NIN,  101 ) 
READ(NIN,  101 ) 


NOROL, NCQUT , NIOUT 

DELI,  DEL2,  TSTOP. ALTST, DELPT 

SGUID, SFRLN 

XIO, YIO, ALTO 

UO,  VO,  WO 

THETO,  PSIO, PHIO 

PQ,  QO,  RQ 

CLP.  CMQ,  CLDA,  DELTF,  8,  SREF 
CLFIN-  DCFIN, CPFIN, PHFIN  FIFAC 
DXCP, FXCN, FXCD 


READ (NIN.  100)  NMACH 
DO  1  I  =  1,  NMACH 

1  READ(NIN,  101)  AMACH ( I ) , CDOT ( I ) , CNAT f I ) ,  XCPT ( I ) 


READ (NIN, 101)  WTAPE, XCGTP, RCGTP, PCGTP.  XIXTP, XIYTP- DIYZ 
READ (NIN, 101)  WPLD, XCGPL- RCGPL. PCGPL 


READ (NIN, 100)  MTYPE 
DO  2  1=1, MTYPE 

READ (NIN,  101 )  WTE ( I ) , WTPR  < I > ,  SP I MP ( I ) ,  CTHR ( I ),  CTHRT (I ' 
READ(NIN, 100)  NTTHR ( I > 

U=NTTHR ( I ) 

DO  3  K  =  i , J 

READ  (NIN,  101  >  TTHRT (K,  I ) .  TTHR (K.  I , 

CONTINUE 


READ(NIN,  100)  NKICK 
DO  A  I “1, NKICK 

READ (NIN,  100)  NTYPE ( I )  ,  NMQT ( I ) 

READ (NIN,  101 )  T I G ( I ) ,  X THR l I ) ,  RTHR  >  I ; ,  PTHR ( I ) ,  DTHR ( I i ,  ATHR <  1 
READ  (NIN,  101  )  XEMPT(  I ) ,  XPROP  ( I ) .  REMPK  I  » ,  RPROP  (  I  ),  PEMPT  I> 
♦PPROP ( I ) 

4  READ (NIN, 101 )  TORQL( I ) , 1CCD ( I ) , TCCL ( I ) , TCCP ( I ) 


READ (NIN, 100)  NALTW 
IF  <  NALTW )  3.3,6 

6  DO  7  1=1, NALTW 

7  READ  (NIN,  101)  ALTWt I > , WIXA< I > • WIYA< I) .  WI2A( I > 


3  REAO(NXN, 100)  NALTA 
IF \ NALTA)  3,  S,  9 


UNCLASSIFIED 


iO  fO 


UNCLASSIFIED 


SM  1081 


9  DO  10  1  =  1,  MALTA 

10  READ  (NIN,101)  ALTA(  I  ) ,  TliPR  <  I  >  #  PRES(  1  ) 

8  CONTINUE 

READtNIN,  100)  NFIN 
IF(NIOUT)  99,99,11 

11  WRITElNOUT, 200) 

WRITE (NCJUT ,  201  )  NOROL,  NCOUT,  NIOUT 

WRITE  (NOUT,  202)  DELI,  DEL2-  TSTOP,  ALTST.  DELPT 

WRITE(NOUT, 203)  SGUID, SFRLN 

WRITE<NQUT,  204)  XIO, YIQ, ALTO 

WRITE(NOUT,  205)  UO, VO, WO 

WRITE(NOUT, 206)  THETO, PSIQ,  PHIO 

WRITE<NOUT,  239)  PO,  00,  RG 

WRITE(NOUT, 207)  CLP,  CMQ, CLDA, DELTF , B, SREF 
WRITE ( NOUT , 221 )  CLFIN,  DDF  IN,  CPFIN,  PHFIN,  FI FAC 
WRITE ( NOUT , 222)  DXCP, FXCN. FXCD 

WR I TE ( NOUT, 208 )  NMACH 
DO  12  1*1, NMACH 

12  WRITE  (NOUT,  209)  AMACH «.  I  > ,  CDOT (  I ) ,  CNAT <  I ) ,  XCPT <  I  ) 

WR I TE ( NOUT ,  210)  WTAPE,  XCSTP, RCGTP, PCGTP, XIXTP.  XIYTP.  DIYZ 
*,  WPLD,  XCGPL,  RCGPL,  PC.GPL 

WRITE (NOUT, 211 )  MTVPE 
DO  13  1  =  1 ,  MTYPE 

WRITE (NOUT, 212)  WTE ( I )  -  WTPR ( I ) , SP IMP < I ) , NTTHR ( I ) , CTHRC I ) , CTHRT ( I ) 
J=NTTHR ( I ) 

DO  14  K=t-  J 

14  WRITE(NOur,  213)  TTHR  T  «  K,  I  ) .  TTHR ( K,  I  > 

13  CONTINUE 

WRITEtNOUT, 214)  NKICK 
DO  13  1=1, NKICK 

15  WRITE  (NOUT,  213)  NTYPEC  I  )  ,  NMOT  (  I  ) ,  TIQ(  I  ),  XTHR  (  I  )  ,  RTHR  (  I  ) ,  PTHR  I  5  , 
*DTHR ( I ) ,  ATHR ( I ) , 

EXEMPT  (  I  ) ,  REMPT  (  I ) .  PEMPT-  I  >,  XPROP I ) ,  RPROP  ( 1 ),  PPROP  <  I  ' 

*TORQL(I>,  TCCD(I), TCCL<I>, TCCP(I) 

WRITEtNOUT, 216)  NALTW 
IF  ( NALTW )  16,  16-  17 
1^  DO  18  1=1, NALTW 

18  WRITE  (NOUT, 217)  ALTW  (  I  ,  WI  XA  (  l  J .  WI YA  (  l>  -  W.IZA  (  I  ) 

16  WRITE  (NOUT, 218)  MALTA 
IF(NALTA)  19,19,20 

C  DO  21  1=1, MALTA 

1  WRITE* NOUT ,  3 1 9 >  a«.TAi  I  • ,  TMPR  (  I  )  ■  PRES «  I  ' 

1 9  CONTINUE 
WRITEtNOUT, 220)  NFIN 

99  CONTINUE 
RETURN 


100  FORMAT <712 ) 

101  FORMAT ( 8F10.  2) 

102  FORMAT (2013) 

103  FORMAT ( 18A4 ) 


UNCLASSIFIED 


104  F  QRM AT ( 3 A4 .  / , 3A4 ,  / , 3A4  > 


UNCLASSIFIED 


SM  1081 


200  FORMAT  <  30 X  ,  '  INPU  T  DATA  " ,  / ,  X ,  1 00  (  IW  1  <  /  ) 

201  FORMAT  IX,  '  NOROL,  NCOUT,  NIOUT  • ,  31 5  ) 

202  FORMAT * X .  'DELI . DEL2, TSTOP, ALTST, DELFT •. 2F10  7, 4F10.  2 ) 

203  FORMAT <  X,  'SGUID, SFRLN  3F10  2> 

204  FORMAT t  X i  ‘  X 10.  YIO- ALTO  '.3F10.2/ 

205  FORMAT  *  X  ,  'UO.  VO. WO  ',3F10.  2> 

206  FORMAT* X, 'THETO, PSIQ, PHIO  ',3F10.2> 

207  FORMAT *  X ,  ' CLP,  CMG,  CLDA,  DELTF,  B,  SREF  ',  6F10.  2.  /) 

221  FORMAT  (  X,  '  CLFIN,  DDF  IN,  CPF  IN,  PHFIN,  FIF  ' ,  5F10.  2,  •' ) 

222  FORMAT  (X,  'DXCP, FXCN, FXCD  ',  3F10.  2.  /  > 

208  FORMAT ix, 'AERODYNAMICS  NMACH= ' , 12, 

*•/,  '  AMACH  CDOT  CNAT  XCPT  ,/> 

209  FORMAT (  X,  F9.  2,  3F10.  4) 

210  FORMAT* /, X, 'WTAPE- XCGTF, RCGTP, PCGTP, X IXTP, XIYTP, DIYZ- ' 

*7F10.  2.  /.  X.  'WPLD,  XCGPL,  RCGPL.  PCGPL  = 

*4F10.  2.  /) 

211  FORMAT ( X ,  'NO.  OF  MOTOR-TYPES  MTYPE-',I2; 

212  FORMAT;/, X,  'WTE. WTPR, SPIMP  , 3F1G  2,  / 

*  X,  'NTHR,  CTHR,  CTHRT  ,110,  2F10.2. 

*/,  X,  'THRUST', 

*/.  '  TTHRT  THRT  '  > 

21 3  FORMAT *F10  4, F10  I) 

214  FORMAT*/, X, 'NO  OF  KICKS  NKICK='. i2, 

*/, '  NTYPE  NMOT  TIG  XTHR  RTHR  PTHR  DTMR  AT HR  XEMF 

*  REMPT  PEMPT  XPROP  RPROP  PPROP  TORGL  TCCD  TCCL  TCCP  / 

215  FORMAT  <215,  3X,  12F7.  2,  F7  5,  F7.  4,  3F7  2) 

216  FORMAT*/, X, 'WIND  NALTW= ' , 13, 

*/,  '  ALTW  WIXA  WIYA  WIZA',/> 

217  FORMAT  <F10.  0*  3F10.  1  ) 

216  FORMAT*/,  X,  'ATMOSPHERIC  DATA  MALTA-',  13, 

*/,  '  ALTA  TEMP  PRES',/) 

219  FORMAT < F 10  0. F10  2,  F10.  3) 

220  FORMAT*/,  X,  'NO.  OF  SHOTS  NFIN«= ' ,  13) 

239  FORMAT*  X,  'PO,  GO,  RO  ',3F10.  2,  /) 

C 

C 

c 

END 


UNCLASSIFIED 


UNCLASSIFIED  SM  1081 

SUBROUTINE  RUNK1 (KUTTA, DT, X, XD, I) 

DIMENSION  C 1  (25 )  ,  C2<  25 .  C3  < 25)  ,  C4 1 25  > ,  SX  ( 25 ) 

COMMON/DEG/C1.  C2,  C3,  C4,  SX 
GO  TO  < 1 . S. 3. 4  >  * KUTTA 

1  SX(I)=X 
C1(I)=XD*DT 
X=SXU)+0.  5*C  1(1) 

RETURN 

2  C2<I>=XD*DT 

X  =  SX  < I )+0. 5*C2( I ) 

RETURN 

3  C3( I )  »  XD*DT 

X  *  SX  < I )+C3( I > 

RETURN 

4  C4<I)  =  XD*DT 

X  =  SX<I )+(Cl ( I >+C4( I )+2. *<C2( I )+C3( I ) ) )/6. 

RETURN 

END 


UNCLASSIFIED 


UNCLASSIFIED  SM  1081 

SUBROUTINE  ATMOS ( RHO,  PP,  TK,  VS,  G,  ALT,  TT,  TP,  TA,  N,  RHFCT,  VSFCT  > 
DIMENSION  TT(IOO), TP<100>, TD(100>, TA(lOO) 

DATA  RE  /20855531  / 

R  *  RE/ <  ALT+RE ) 

H  =  .  3048*ALT*R 

TK  «  288.  13-0.  0065*H 

PN  =  101300.  *<TK/288  1 3  >  **  t -*-5  255  > 

IF <  N.  LE  1 )  GOTO  3 

CALL  INTEGITT,  TP,  TD-  TA,  DT,  DP,  D,  ALT,  N) 

TK  =  TK#<  1.  +DT/100.  ) 

PN  =  PN*<  1.  +DT/100.  ) 

3  CONTINUE 

PP  a  0  G20S85*PN 
RHO  »  6  75976E-6*PN/TK 
VS  =  65.  7688*SQRT<TK) 

G  *  32.  1741*R**2 
RHO  »  RHO*RHFCT 
VS  =  VS* VSFCT 
RETURN 
END 


UNCLASSIFIED 


UNCLASSIFIED  SM  1081 

SUBROUTINE  INTEQ  <TY1>  TY2,  TY3i  TX,  Yl,  Y2,  Y3,  X.  N) 

Nl-N-1 

IF < X.  LT.  TX<  1 )  )  GO  TO  3 
IF< X.  GT.  TX (N)  )  GO  TO  4 
DO  1  1*1,  N1 
J=I 

IF ( ( X-TX ( I ) )*( X-TX (1+1)))  2,2,1 

CONTINUE 

CONTINUE 

J*N1 

GO  TO  2 

J*1 

CONTINUE 

DTY1*TY1 < J+l )-TYl ( J) 

DT Y2*T Y2  <  J+ 1 ) -T Y2 ( J ) 

DTY3=TY3 <  J+l  > -TY3  <  J ) 

DX*(X-TX< J) ) / ( TX ( J+l ) -TX ( J) ) 

Y1=TY1 < J)+DTY1*DX 
Y2«=TY2  ( J )  +DTY2*DX 
Y3*TY3  <  J  >  +DT Y3*D  X 
RETURN 
END 


UNCLASSIFIED 


UNCLASSIFIED  SM  1081 

SUBROUTINE  UlINDtTX,  TV,  TZ,  TA,  N,  ALT,  A,  Rl.  R2,  R3  ’ 

DIMENSION  TX(N),  TY(N),  TZ(N),  TA<N),  A <3.  3),  R(3>.  X(3) 

CALL  INTEQ(TX,  TY,  TZ.  TA,  X  (  1 ) ,  X  (2 > ,  X  <3 )  <  ALT,  N> 

DO  1  1  =  1,3 
R ( I ) =0 
DO  1  J-- 1 , 3 

R  < I ) *R  < I >  +  X< J)*Al I, U> 

CONTINUE 
R1*R<  t  > 

R2=R ( 2  > 

R3=R ( 3 > 

RETURN 

END 


UNCLASSIFIED 


UNCLASSIFIED  SM  1081 

SUBROUTINE  INTPL  < TY.  TX,  Y,  X,  K.  N > 

DIMENSION  TX(N),TY<N> 

N1  =  N-  1 

DO  1  I  =  K,  N1 
J=I 

IF< ( X-TX ( I ) )*  < X-TX ( 1+1 > ) ~0  OOOl  >  2.2,1 

1  CONTINUE 
U=N1 

2  CONTINUE 

DTY  =  TYlU+l)  -  TY ( J ) 

Y  =  TY ( J )  +  DTY* (X-TX(U) ) / ( TX ( J+l >  — TX ( J ) ) 

K— J 

RETURN 

END 


UNCLASSIFIED 


o  n  r>  o  o  n  n 


UNCLASSIFIED  SM  1081 

SUBROUTINE  OUTPU  (A. IQUT, NQUTi NCQUT, NSHOT ) 

C 

DIMENSION  A(60, 200) 

DATA  Si  P,  F  /l.  354,  .  453,  .  3048/ 

NC0NV=NC0UT/10 
NTABL=NC0UT-NC0NV*10 
IF  <NTABL-1 )  99,  1,  1 
1  CONTINUE 

1ST  TABLE 

WRITE ( NOUT , 100)  NSHOT 
IF(NCONV)  3,3,2 

3  WRITE (NOUT , 101 ? 

GOTO  4 

£  WRITE (NOUT,  102) 

DO  8  1=1, I OUT 
DO  7  J=2,  5 

7  A  (  J,  I  )  =A  (  J,  I  )  *F 
DO  6  J=S, 10 

6  A<  J,  I  )=A<  J.  I  >*F 

8  CONTINUE 

4  CONTINUE 
DO  5  1=1, [OUT 

5  WR  ITE  ( NOUT ,  103)  ( A  <  J,  I  )-  J=l,  8),  A-I59,  I  ),  A(6Q.  I  ) 

IF  (NTABL-2)  99,  10,  10 
2ND  TABLE 

10  WRITE  (NOUT, 100)  NSHOT 
WRITE (NOUT  >  110> 

DO  1 1  1=1, I OUT 

11  WR  ITE  (NOUT,  111  )  A(l,  I),  <A(J,  I),U=11,  19) 

/= 

S_r 

IF  (  NTA6L-3  )  99,20.20 

C 

C  3RD  TABLE 


r 

C 

r 


20  WRITE (NOUT, 100)  NSHOT 
WRITE (NOUT, 120) 

IF(NCONV)  21,21,22 

21  WRITE(NOuf,  121) 

GOTO  24 

22  WRITE ( NOU T ,  122) 

DO  27  1=1, I OUT 

A (22,  I ) =A ( 22,  I )*P*9.  8065/F**2 
A (36,  I  )  =A ( 36,  I  )*P 
A  (37.  I  >-=A(37,  I  )*F 
A ( 38,  I >  =A ( 33,  I  )*F 
DO  27  U=53, 55 

27  A(  J,  I  )  =  A(U,  I  >*F 
24  DO  28  1=1, IOUT 

28  WRITER  NOUT,  123)  A  (  1 ,  I  ) ,  (  A  t  J,  1  ) ,  J=20,  23  )  ,  <  A  (  J,  I  ,  J  --36  36  > 
*  (  A  (  J,  I  )  .  J=53,  55) 


IF ( NTABL-4 )  99, 30,  30 
4TH  TABLE 

■ 

aSHOT 


30  WRITE (NOUT, 100) 


UNCLASSIFIED 


UNCLASSIFIED 
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UR I TE ( NOUT > 130) 

IF(NCONV)  31 <  31.  32 

31  UR  I TE ( NOUT .131) 

GOTO  36 

32  WRITE (NOUT,  132) 

DO  33  1*1, I OUT 
DO  34  U=24,  35 

34  A(J,  I)=A(J,  I)*P*9.  8065 
DO  35  J=30,  35 

35  A  ( J,  I  )=A(  J,  I  )*F 

33  CONTINUE 

36  DO  37  1=1,  I OUT 

37  WRITE  (NOUT,  133)  A(l,  I),  ( A  ( J,  I),  J=24,  35 ) 

IF(NTABL-S)  99,40,40 

5TH  TABLE 

40  WRITE (NOUT, 100)  NSHOT 
WRITE (NOUT,  140) 

IF (NCONV )  41,41,42 

41  WR I TE ( NOUT  ,141) 

GOTO  47 

42  WR I TE ( NOUT ,  142) 

DO  43  1  =  1,  I OUT 
DO  45  J  =  56, 58 

45  A( J,  I)  =  A( J,  I)  *  F 
DO  46  J=47, 52 

46  A(J, I)  =  A ( J, I )  *  P*F*9  8065 

43  CONTINUE 

47  DO  48  I  =  1,  I  OUT 

48  WRITE (NOUT,  143)  A< 1,  I ) ,  ( A< J.  I > , J»47. 52) ,  ( At J,  I > , 0=36,  58) 
IF  ( NTABL-6 )  99,50,50 

6TH  TABLE 

50  WRITE (NOUT,  100)  NSHOT 
WRITE (NOUT,  150) 

IF  (NCONV)  51,51,52 

51  WR ITE( NOUT, 151 ) 

GOTO  58 

52  WRITEtNOUT,  152) 

DO  53  1  =  1,  I OUT 

A(41, I)  =  A(41, I)  #  F 
A(42, I)  =  A ( 42, I)  *  F 
A (39,  I  )  =  A(39,  I  >  *  S 
A  (40,  I  )  =  A  (40,  I )  *  S 
DO  54  U=43, 46 
54  A(  J.  I )  =  A(  J,  I )  *  S 

53  CONTINUE 

58  DO  59  1  =  1,  I OUT 

59  WRITE  (NOUT,  153)  A(l,  I),  ( A  (  J,  I),  J=36,  37).  A(41,  I ),  A  (42,  I ) , 
*A(40,  I),A(39,  I),  (A(J,  I  ) ,  J=43,  46) 


99  CONTINUE 
RETURN 

100  FORMAT <1H1,/,X,40(1H*>,  '  RESULTS  SHOT  NO.  ',  12,  X,  40(  1H*> , // ) 

101  FORMAT  ( '  TIME  RANGE  CRpSS-R  ALT  VELOC 

UNCLASSIFIED 
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*V  AZIM  U-DOT  ASP -EL  V  ASP -A  l',/, 

*  •  (SEC)  ( FT )  (FT)  (FT)  (FT/S>  ( 

*G>  (DEG)  (FT/S2)  (DEG)  (DEG)  ',/> 

102  FORMAT  f  TIME  RANGE  CROSS-R  ALT  VElOC  E 

*V  AZIM  U-DOT  ASP-ELEV  ASP-AZ',/, 

(SEC)  ( M )  < M )  <M)  (M/S)  ( 

*G>  (DEG)  (M/32)  (DEG)  (DEG)'. /> 

102  FORMAT  <F10.  2.  3F10.  0,  6FlQ  2! 

110  FORMAT  ('  TIME  THETA  'PSI  PHI  P 

*  R  P-DOT  G-DOT  RDOT ' ,  /, 

*  '  (SEC)  (DEG)  (DEG!  (DEG'-  (DEG/S)  (DE 

•>L.  (DEG/S)  (DEG/32)  (DEG/S2)  (DES/S2J',/) 

111  FORMAT  1 OF  10.  2) 

120  FORMAT  ; '  TIME  ALPHA  BETA  DYN. PRESS  MACH  W 

*T  C  G  C.  P  WIND-X  WIND-Y  WIND-Z  '  ? 

121  FORMAT  ('  (SEC)  (DEG)  (DEG)  (LB/FT2)  (-)  ( 

•*)  (FT)  (FT)  (FT/S?  CFT/S)  (FT/S)'./) 

122  FORMAT  ■  '  (SEC)  (DEG)  (DEG)  (N/M2)  (-.'  ( 

*■>  ( M  >  (M)  (M/S)  (M/S)  (M/S)'./) 

123  FORMAT  (11F10. 2) 

130  FORMAT) •  TIME', 

*  '  FX-AERO  FV-AERO  FZ-AERu  FX-THR  FY-THR  FZ-THR  TX 

*ERO  TY-AERO  TZ-AERO  TX-THR  TY-THR  TZ-THR ' ) 

131  FORMAT ' '  (SEC)', 

*  '  (LB>  (LB)  (LB)  (LB)  (LB)  (LB)  (L 

*FT>  (LB-fi-FT)  (L8-*FT )  <LB*FT>  (LB*FT)  ( LB*FT )  ' ,  /  ) 

132  FORMAT ( '  (SEC) ', 


*  '  ( N )  ( N )  ( N )  ( N )  ( N ) 

*M )  <  N-a-M  )  ( N*M  >  l  N«M  >  ( N*M  >  ( N*M  )  ' ,  /  ) 

133  FORMAT  (F9  2,  F9.  0,  2F9.  2,  F6  0,  2FS.  1, 3F9  1,  3F8.  1  ) 

( N ) 

140  FORMAT ( 

TIME 

COUP-L 

COUP-M 

COUP-N 

DAMP-L 

*  '  DAMP-M 

DAMP-N 

WW1 

WW2 

WW3 

) 

141  FORMAT ( 

(SEC  ) 

(LB*FT) 

(LB*FT> 

(LB*FT) 

(LB*FT> 

*  '  <  LB*FT ) 

(LB*FT) 

(FPS) 

(FPS) 

(FPS) 

,  /  ) 

142  FORMAT ( 

(SEC) 

(N*M> 

(N*M) 

(N*M, 

(  N*M ) 

■a  '  <N*M) 

143  FORMAT ( 1QF10. 2) 

( N*M ) 

(M/S) 

(M/S) 

(M/S) 

,  /  ) 

1 50  FORMAT ( 

TIME 

WEIGHT 

X-CG 

Y-CG 

Z-CG 

*  '  I-XX 

I-YY 

I-ZZ 

I  —  XY 

I-XZ 

I-YZ  ' 

131  FORMAT ( 

(SEC  ) 

(LB  ) 

(FT) 

(FT) 

(FT) 

*  '  (SL*FT2> 

■*/  > 

(SL*FT2) 

(SL*FT2> ( SL*FT2 ) 

(SL*FT2> 

<SL*FT2) ' 

152  FORMAT ( 

(SEC  ) 

(KG) 

(M) 

<  M ) 

(M) 

*  '  <KG*M2> 

(KG*M2) 

(KG*M2) 

(KG*M2> 

(KG*M2) 

<KG*M2> 

*) 

153  FORMAT ( 1 1F10.  2) 


END 


UNCLASSIFIED 


UNCLASSIFIED  SM  1081 

SUBROUTINE  OUTAB  (NOUT,  NCOUT,  NSHOT.  RANGE.  FLTIM.  APOGE,  TAPOO.  QBMAX. 
♦XMMAX.K) 

DIMENSION  RANGE ( 100),  FLTIM! 100).  APOQE! 100), TAPOO! 100),  OBMAX(IOO), 
♦XMMAX  < 100) 

DATA  F, P  /.  3048,  .  453/ 

GOTO  < 1, 2) >  K 

1  I “NSHOT 

WRITE (NOUT, 100)  I, RANGE! I ) , FLTIM ( I ), APOOE! I). TAPOO! I),  XMMAX!I). 
♦QBMAX ! I > 

RETURN 

2  WRITE ! NOUT , 110) 

IF !NC OUT-10)  3,3,4 

4  WRITE! NOUT,  114) 

DO  3  1*1, NSHOT 
RANGE! I)  *  RANGE ! I ) #F 
APOGE  < I )  *  APOGE ! I ) #F 

5  QBMAX!  I)  =  QBMAX!  I  >♦?♦<?.  8065/F##2 
GOTO  6 

3  WRITE! NOUT,  115) 

6  DO  7  1*1, NSHOT 

7  WRITE!NOUT,  120)  I , RANGE ! I > , APOGE ! I > , FLTIM! I ) , TAPOO! I ) , XMMAX ! I > , 
♦QBMAX! I > 

RETURN 


100  FORMAT!///,  X.  10!1H#>,  '  SUMMARY: NO, RANGE, FLICHTTIME. APOGEE.  T-APOOEE 
*,  MAX.  MACH,  MAX.  DYN.  PRESS'  .  /,  122.  FB  0,  F8.  1.  F8.  O.  F8.  1,  F8.  2.  F14.  1  ) 

110  FORMAT < 1 HI ,  X,  40 ( 1H* ) .  '  SUMMARY-TABLE  ' , 40! 1H^ > . / / ,  X, 

♦  'NO.  ' ,  5X,  'RANGE'.  4X.  'APOGEE',  X,  'FLIGHT-TIME '.  X,  'T-APOOEE ' .  X. 

♦  'MAX.  MACH'.  X,  'MAX.  DYN.  PRES') 

114  FORMAT!  10X,  '(M)  '.  6X,  '!M>  ',  6X,  '!SEC)  ',  5X.  '(SEC)  5X,  ' !  1  '.  7X.  '(N/M2 

♦  )',/) 

115  FORMAT!  lOX.  '(FT)  ',  5X.  '  (FT)  ',  5X,  '  (SEC)  '.  5X,  '  (SEC)  ',  5X,  '  (  >  ',  6X.  '(LB 
*/FT2 )  ',  /) 

120  FORMAT!  13,  Fll.  0,  F9.  0,  Fll.  2,  F10.  2,  F7.  2,  F13.  1) 

END 
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APPENDIX  2 


THRUST  FORCE  AND  MOMENTS  RESOLVED  IN  F* 


In  Section  2.6  of  the  main  text,  the  thrust  forces  and  moments  were  presented  in 
terms  of  components  in  the  body- fixed  reference  frame  F„.  In  order  to  obtain  these 
values,  however,  a  transformation  is  required  taking  the  thrust  vector  and  transforming 
it  into  components  in  Fa  and  then  using  these  components  and  the  location  of  the  point 
of  action  of  the  thrust  vector  relative  to  F#,  to  compute  the  thrust  moments.  This 
Appendix  summarizes  this  analysis. 

Only  the  case  for  one  rocket  motor  of  arbitrary  orientation  relative  to  the  fuselage 
reference  line  (FRL)  is  considered.  The  results  for  several  rocket  motors  follow  with  an 
algebraic  summation  of  each  motor’s  contributions. 

A  reference  frame  Fr##,  is  defined  in  Figure  A2. 1 .  Frw.  has  its  origin  located  at  the 
point  of  action  on  the  airframe  of  the  thrust  vector  of  the  i-th  rocket  motor,  and  its 
x-axis  in  the  negative  thrust  direction.  The  orientation  of  the  axes  of  FTh,  relative  to  axes 
of  the  structural  reference  frame  F*  is  given  by  the  ordered  Euler  angle  rotations  (4™,, 
Qth,,  dr##,),  i.e.  a  rotation  $rw<  about  the  x-axis  of  F«,  a  rotation  BTHi  about  the  y-axis  of 
the  reference  frame  resulting  from  the  first  rotation,  and  a  rotation  6th,  about  the  z-axis 
of  the  reference  resulting  from  the  second  rotation.  A  further,  important  constraint  is 
that  the  rotation  $THi  corresponds  to  the  angular  cylindrical  coordinate  of  the  point  of 
action  of  the  thrust  vector  of  the  i-th  rocket  motor  relative  to  the  FRL  (see  Figures  4  and 
A2.1). 


In  terms  of  the  data  inputs  to  the  BALS1M  package,  the  Euler  angles  0TH„ 
6th,)  have  a  one-to-one  correspondence  with  the  program  variables  (PTHR(I), 
DTHR(I),  ATHR(I))  (see  Appendix  1). 

Since  the  orientation  of  Ffl  is  related  to  the  orientation  of  F*  by  a  n  rotation  about 
the  y-axis  of  F«,  it  may  be  shown  that  the  rotation  matrix  L BTHi  relating  components  of  a 
vector  expressed  in  FTH<  to  the  components  of  the  same  vector  expressed  in  F„  is  given  by: 
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—  COSdrH,COS0rH, 

sindrW(cos0ra< 

-  sin0r«, 

L  BTHj  ~ 

cos  d  th,  sin0rH,  sin^TH,. 

+  Sindr//,COS^r«j 

-  sindr^sinSrHiSin^™, 

+  COS  d  th,  COS  ^  th, 

—  cos  0rH,sin<|>  TW< 

(A2,l) 

cos  d  THt  si  n  0  th,  cos  $  th, 

—  sindrz/jSin^rHj 

—  sin  d  th,  sin  0  th  ,  cos  ^  th, 
cosd™(sin<|>TWi 

—  COS  0  th,  COS  ^  th, 

From  the  definition  of  Fth,  it  follows  that  the  thrust  of  the  i-th  rocket  motor  is 
given  by: 


T  ™, 

=  (  -T,,  0,  0)r 

(A2,2) 

and  the  nozzle  moment  of  the  i-th  rocket  motor  is  given  by: 

=  ( —  Lnij,  0,  0)r 

(A2,3) 

Using  (A2,l)  and  summing  over  all  the  rocket  motors,  it  follows  that: 

Xtb 

Mm 

=  ,  Z  T,  COS  drn,  COS  0th, 

1=1 

(A2,4a) 

Yr, 

Nm 

=  .  -T,(cosdTH,sin0TH,sin«J>TH,  +  sindrH,cos<|>rH,) 

(A2,4b) 

Zr a 

Nm 

=  .  2 1  — T,(cosdrH, sin0rH( cos^th,  —  sindm, sin^rn,) 

(A2,4c) 

and 

(Ltb)„ 

Mm 

=  j^i  COS0T-/T, 

(A2,5a) 

(Mt,)„ 

Nm 

=  .  Z) -L„,(cosdrH,sin0TH,sin<j>TH,  +  sindrH,cos«|>rH,) 

(A2,5b) 
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(NTb)„  =  Z  -L^CcosdTH.sinerH.cos^rH,  -  sind™, sin^ra,)  (A2,5c) 

What  remains  is  to  compute  the  thrust  moment  (h£n)««  generated  by  the  i-th 
rocket  motor  due  to  the  displacement  of  its  point  of  action  relative  to  the  vehicle  centre- 
of-mass.  Let  r^TH.  be  the  vector  from  the  vehicle  centre- of- mass  to  the  origin  of  Frw,  (see 
Figure  A2.1),  (xr«„  y Th„  z th.)  are  the  coordinates  of  the  origin  of  FTHi,  in  F*,  and  (xc„ 
ycg,  zcg)  are  the  coordinates  of  the  vehicle  centre- of- mass  in  F«.  *s  given  by; 

(Mr,)c*  =  I+THi  x  (A 2,6) 

or 

Mr,  =  1*1?  (A2,7) 


where  r  *  is  the  matrix  equivalent  of  the  vector  cross-product  (Reference  3).  Using  the 
definitions  of  the  previous  paragraph  and  Figure  A2.1,  it  may  be  shown  that: 


0 

(Z THt  ~  Z eg) 

(y™,  - 

y  eg) 

-(Zr„.  -  Zct) 

0 

(Xth,  — 

XCg) 

-(y™,  -  y«) 

-(Xrw,  -  xc,) 

0 

From  (A2,4),  (A2,7)  and  (A2,8),  it  follows  that: 


(A2,8) 


N„ 

(Lrs)c*  —  .2  l— (z m  ~  zc<) T , (cosdw, sin@TW. sin^r«,  +  sind7-w<cos$7-j/.) 

-(y THt  -  yJT^cosdr^siner^cos^rH,  -  sindTWi sin<j>TO,)] 


(A2,9a) 


nm 

(Mrs)c*  =  [  (z THt  ~  zc,)T,cosdw,coserH< 


(A2,9b) 


-(Xr  -  xcm)TI(cosdm,sin0rH/cos«j>rH,  -  sindre,  sin^w,)! 
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N* 

(Nrs)c,  =  .  I  [ -(yr//,  -  ycJT^osdr^.cose™,  (A2,9c) 

1=1 

+  (Xr  -  xcJT<(cosdr*,sin8r«,.cos<j>r#,1  +  smdrHi COS <!>„/,)] 

All  the  terms  that  are  required  to  specify  the  total  thrust  moments  (Lra,  MT«,  Nril) 
as  given  by  equations  (2.6,1a),  (2.6,1b)  and  (2.6,1c)  have  now  been  specified. 
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